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In this Chapter we provide an overview of the current first-principles perspective on 
flexoelectric effects in crystalline solids. We base our theoretical formalism on the 
long-wave expansion of the electrical response of a crystal to an acoustic phonon 
perturbation. In particular, we recover the known expression for the piezoelectric 
tensor from the response at first order in wavevector q, and then obtain the flex¬ 
oelectric tensor by extending the formalism to second order in q. We put special 
emphasis on the issue of surface effects, which we first analyze heuristically, and 
then treat more carefully by presenting a general theory of the microscopic response 
to an arbitrary inhomogeneous strain. We demonstrate our approach by present¬ 
ing a full calculation of the flexoelectric response of a SrTiOs film, where we point 
out an unusually strong dependence of the bending-induced open-circuit voltage on 
the choice of surface termination. Finally, we briefly discuss some remaining open 
issues concerning the methodology and some promising areas for future research. 


1. Introduction 

First-principles electronic structure calculations have played an increasingly impor¬ 
tant role in our understanding of the properties of materials and nanostructures in 
recent decades. The phrase “first principles” is generally used in the condensed-matter 
community to convey the notion that the calculations are free of adjustable parame¬ 
ters, taking as input only some list of atoms, their atomic numbers, and some initial 
guesses at their coordinates in the unit cell. One then solves the Schrodinger equation 
for the electrons in some approximation, computes the relaxed atomic coordinates, 
and calculates the desired properties of the crystal. In the condensed-matter commu¬ 
nity this is typically done in the framework of density-functional theory (DFT),li]as 
shall be assumed below, but Hartree-Fock or other quantum-chemical methods can 
also be used. 
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While the accuracy and efficiency of DFT methods have improved over the years, 
of equal importance has been the increasing range of quantities that can be computed. 
In the context of dielectric properties, the implementation of linear-response theory 
for phonon and electric-field perturbations in the 1980s and 1990s opened up the 
calculation of phonon frequencies, dynamical charges, and both electronic and lattice 
contributions to the dielectric constant.^ While there was initially some doubt about 
whether the piezoelectric response was a bulk property at all, a seminal paper of 
Martin laid this question to rest,l^ and the computation of the piezoelectric tensor 
is now a standard feature of most DFT codes as well. Strangely, although many 
of the above properties can be computed as derivatives of the electric polarization 
P, a proper definition of the polarization P itself proved more difficult; the physics 
was clarified, and practical methods for computing it, were developed only in the 
mid-1990s with the appearance of the “modern theory of polarization.”®^ Related 
methods for computing the orbital magnetization of ferromagnets and the properties 
of crystals in finite electric fields have been developed since the 2000’s.l3 

The flexoelectric tensor has been among the few physical properties to have re¬ 
sisted a proper first-principles formulation even until today. The theory of flexoelec- 
tricity was pioneered in the 1980s by TagantsevEEl However, because it encodes a 
response to a strain gradient, rather than just a strain, and because a strain gradi¬ 
ent is inconsistent with ordinary cell-periodic boundary conditions, methods based on 
Bloch’s theorem cannot be straightforwardly applied in the first-principles context. A 
serious attack on this problem did not begin until 2010, when Hong and collaborators 
presented the results of calculations on supercell configurations containing strain gra¬ 
dients.^ Subsequent papers of ResteP^^ and Hong and Vanderbilt?^ clarified aspects 
of the electronic contribution to the flexoelectric response. 

More recently, StengeP^ and Hong and VanderbiltP^ tackled the problem in a sys¬ 
tematic way, and working from slightly different perspectives, arrived consistently at 
a nearly complete framework for defining, and eventually computing, the flexoelec¬ 
tric tensors fully from first-principles. Some components of the flexoelectric tensor 
that can be expressed only in terms of bulk current responses, as opposed to charge 
responses, still require care in their interpretation and await the development of ef¬ 
ficient methods for calculating them. However, we can expect these difficulties to 
be cleared up soon, so we can look forward to a new era in which first-principles 
calculations of flexoelectric responses can flourish and contribute to a fast-evolving 
experimental field. The purpose of this chapter is to outline the physical principles 
underlying these advances in the understanding and computation of flexoelectric re¬ 
sponses, and to summarize a few of the preliminary results that have been presented 
in the literature to date. 

2. Theory and methods 

2.1. Strain, strain gradients, and responses 

We begin by establishing our notation. In continuum mechanics, a deformation can be 
expressed as a three-dimensional (3D) vector field, UQ,(r), describing the displacement 
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of a material point from its reference position at r to its current location r',|^ 

r^(r) =ra+ Ua{r), 


The deformation gradient is defined as the gradient of Ua taken in the reference 
configuration, 

- , , . , dua{r) 

eaft{v) = Ua,0{v)= . ( 1 ) 

ectp{r) is often indicated in the literature as unsymmetrized strain tensor”, as it 
generally contains a proper strain plus a rotation. By symmetrizing its indices one 
can remove the rotational component, thus obtaining the symmetrized strain tensor 

^a/3 — 2 • 


This eafi is a convenient measure of local strain, as it only depends on relative dis¬ 
placements of two adjacent material points, and not on their absolute translation or 
rotation with respect to some reference configuration. 

In this work we shall be primarily concerned with the effects of a spatially inhomo¬ 
geneous strain. The third-rank strain gradient tensor can be defined in two different 
ways, both important for the derivations that follow. The first (type-I) form consists 
in the gradient of the unsymmetrized strain. 




dSapir) 

dr^ 


d‘^Ua{v) 

drpdr^ 


( 2 ) 


Note that rja^g-y, manifestly invariant upon ft j exchange, corresponds to the Vafij 
tensor of Ref. 12 and to the symbol deap/dr^ of Ref. Alternatively, the strain 
gradient tensor can be defined (type-II) as the gradient of the symmetric strain, Eq,^, 




dsap{r) 

dr^ 


invariant upon a ft exchange. It is straightforward to verify that the two tensors 
contain exactly the same number of independent entries, and that a one-to-one rela¬ 
tionship can be established to express the former as a function of the latter and vice 
versa. 


Va,P'y ^a/3,7 


^Py,a- 


(3) 


The piezoelectric and flexoelectric tensors describe, respectively, the macroscopic 
polarization response to a uniform strain and to a strain gradient. In type-I form, 
these are 

_ dPa 

6a/37 — 7 5 

I _ 


(4) 

(5) 


^As in this Chapter we shall deal exclusively with linear flexoelectricity, we shall assume a regime 
of small deformations henceforth. 
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While the type-I form is more convenient to derive and calculate, the type-II repre¬ 
sentation is often preferred in applications. The type-II flexoelectric tensor is defined 
as 


/^aA,/37 


dP„ 


de 




( 6 ) 


Note that and are both symmetric under the last two indices, and are related 
to each other via Eq. (§ according to 

A^aA,/97 Ma7,A/3 (^) 


Mq/3,7A q (/^aA,/37 ^aj,fiiX^ ■ 


( 8 ) 


2.2. Long-wave approach 

A macroscopic strain gradient breaks the translational symmetry of the crystal lattice. 
For this reason, the response to such a perturbation cannot be straightforwardly 
represented in periodic boundary conditions. This makes the theoretical study of 
flexoelectricity more challenging than other forms of electromechanical couplings such 
as piezoelectricity. To circumvent this difficulty, we shall base our analysis on the 
study of long-wavelength acoustic phonons. These perturbations, while generally 
incommensurate with the crystal lattice, can be conveniently described in terms of 
functions that are lattice-periodic, and therefore are formally and computationally 
very advantageous.^ 

Consider a crystal lattice spanned by the real-space translation vectors R; and 
by the basis vectors in such a way that = R; -I- Tk indicates the location of 
the atom of sublattice k and cell 1. In full generality, the atomic displacements along 
the Cartesian direction a associated with a phonon eigenmode of wavevector q can 
be written as 






(9) 


where (independent of either I or t) is an eigenvector of the dynamical matrix at 
q, and uj is the frequency. 

A convenient description of arbitrary mechanical deformations can be established 
by choosing an acoustic phonon branch, and by performing a long-wave (small q) 
expansion of its eigenvector in the vicinity of the T point. Provided that the long- 
range electrostatic fields are adequately screened (see Sec. 2.3 for a discussion), the 
aforementioned expansion can be written as 


^KCt ap-i ^J^X^aP-fX T ■ ■ ■) ) (19) 

where U is a Cartesian vector, 6ap is the Kronecker delta, and T)^^^ and are 

third- and fourth-order tensors, respectively. (The dots stand for higher-order terms, 
which are irrelevant in the context of the phenomena described here.) At order zero 
in q the phonon eigenmode is a rigid translation of the whole lattice along U (note 
the absence of a sublattice index), while the first- and second-order terms describe 
the internal-strain response of the lattice to a uniform strain or to a macroscopic 
strain gradient, respectively. 
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Of course, to obtain the relevant electromechanical coupling coefficients, the sole 
knowledge of the lattice distortions is not sufficient - one needs to establish a link 
between atomic displacements and macroscopic polarization. While in a simplified 
point-charge model such a link is straightforward, in the case of a more realistic 
quantum-mechanical description of a solid things are significantly more involved, as 
one needs to understand how the electronic wavefunctions, and not only the nuclei, 
respond to a macroscopic deformation. If the deformation is sufficiently slow, which 
is the case of the phenomena described in this chapter, the electronic cloud responds 
adiabatically to atomic motion by generating a microscopic current density (i.e., the 
quantum-mechanical probability current). For example, if we displace by hand one 
atomic sublattice as 




( 11 ) 


the microscopic current density that is linearly induced by such a perturbation can 
be written as0 


J(r,t) = A(t)P^^(r)e*'iF (12) 

The function P)?o(r) is the microscopic polarization response; its cell average, 




(13) 


where is the cell volume, describes the contribution of atomic motion to the macro¬ 
scopic polarization, which is the quantity we are ultimately interested in. 

To go from here to the electromechanical tensors we need one more step, i.e., the 
small-q expansion of P„^. Again expanding in powers of q and keeping terms up to 
second order ,[^ 


P** — P^*^^ in 


o 


(14) 


The zero-th order term is the macroscopic polarization response to a macroscopic 
translation of the sublattice k along the direction /3. This corresponds precisely to 
the definition of the Born dynamical charge tensor Z*, 


p(o) ^ 


(15) 


The remaining P-tensors can be regarded as higher-order counterparts of the Born 
charges. (Physically they are directly related to the moments of the current density 
induced by the displacement of an isolated atomii^lHI) 

Multiplying the lattice-polarization coupling tensors with the phonon eigendis- 
placements, we can collect terms order-by-order in q. The zero-order term (rigid 

'^Recall that, in classical electrostatics, the density of bound currents J and the microscopic polar¬ 
ization P are related by J = dF/dt. _ _ 

‘^Note the difference in sign convention between Eq. and Eq. . In the former case, the choice 
of the sign was uniquely determined by the interpretation of T and N as internal-strain response 
tensors. In the latter case, the adopted convention allows one to identify the tensors with the 

real-space moments of the current-density response 1133 
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translation) vanishes due to the acoustic sum rule. At first order in q, we obtain the 
explicit expression for the piezoelectric tensoi^ 




_ V^p(1.7) 


z: 




(16) 


where the first and second terms are the electronic (frozen-ion) and lattice-mediated 
terms respectively-!^ Collecting the terms at second order in q gives the flexoelectric 
response, which is again a sum 


I _ -I . I,mix I I,latt 

/^a/9,7A ' /^a/3,7A ^a^,'yX'> 

of electronic and lattice terms 

I ^ i V 

2 / ^ a,K.^ i 

K 

I,mix 1 

1^01(3,-yX 2 \ p(3-y-^ oi,Kp ^ p/SX-^ a,Kp J ’ 

I,latt _ ^K,ap 
r^Q:/ 3 , 7 A ^ ^^p0yXi 


(17) 

(18) 

(19) 

( 20 ) 


where the bar symbol on the first term indicates a purely electronic response and ‘mix’ 
and ‘latt’ refer to “mixed” and “lattice-mediated” contributions, respectively. While 
the piezoelectric and flexoelectric responses have been developed in parallel until 
now, we will henceforth concentrate on the latter, referring the reader to Refs. 13|14 


for the detailed treatment of the piezoelectric response. The corresponding type-ll 
flexoelectric responses are 


II _ -II I II,mix I II,latt 

Mq:A,/37 l^OiX,(3y ' /^q:A,/37 /^q:A,/37’ 


where 


MaA 




p(2.7A) p(2,A/3) _p(2./37)\ 

;, k /3 ' a.,Ky ^ ol^kX j 


11,1. 

^aA,/37 

II,latt 
A^aA,/37 


p(iA) 

- p^'y^ a.,K,p 1 


z* 




For later convenience we rewrite Eq. (p4|) as 




^_II,latt ^K^Ctp jK 

f-^aX,l3-y ^2 


( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 


where (the type-Il counterpart of the type-1 internal-strain tensor N) is the 

quantity in parentheses on the right-hand side of Eq. (24). 

To summarize, according to Eq. 0 a long-wavelength sound wave is comprised of 
a lattice-periodic distortion pattern modulated by a time- and space-dependent 
complex phase factor. At zero order in q the deformation can be described as purely 
“elastic,” but at higher orders (i.e., when moving away from the zone center), internal 


The unsymmetrized strain is ep^{r) = iUpq^e^'^ ''\ this can be replaced with the symmetrized 
strain tensor after observing that both terms on the right-hand side are invariant with respect to 
fd'y exchange. 
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relaxations of the basis atoms in the primitive cell occur, as described by the tensors 
r and N (or L) at first and second orders in q, respectively. These are related 
to how the crystal locally responds to a macroscopic strain (first order, “piezo”) or 
strain gradient (second order, “flexo”). The reader is referred to Refs. 13|14 for the 
derivation of explicit expressions for these tensors, but we shall highlight the main 


conceptual issues associated with them in Sec. 2.4 Each consecutive order in Eq. (10) 


gives rise to a corresponding term in the expressions for the flexoelectric tensor in 
Eqs. ( [T7| ) and (21). 

Regarding the purely electronic term, /i, which is associated with the purely 
elastic part (order-zero in q, also re ferre d to as “frozen ion deformation”), we de¬ 
fer its detailed discussion to Section 


2.5 


It can be shown that the “mixed” 
term involving is active only in crystals that are characterized by Raman-active 
phonons,^ which is not the case for simple systems such as cubic rocksalt or per- 
ovskite crystals. (Again, we refer the reader to Refs. 13|14 for the explicit discussion 
of this term.) By contrast, the term is present in any insulator with IR-active 

phonons; as this term is very important in practical applications of the flexoelectric 
effect, we shall discuss it shortly in Sec. |2.4[ First, however, we shall briefly com¬ 
ment on an important issue that is relevant to the above discussion, concerning the 
treatment of the macroscopic electric fields in the long-wave phonon analysis. 


2.3. Macroscopic electric fields 

Depending on their polarity, long-wave phonons in a crystalline insulator generally 
produce macroscopic electric fields. These are due to the charge perturbation that is 
generated by the lattice distortion, and have a nonanalytic behavior in the vicinity of 


the r point. For example, for a monochromatic perturbation such as that of Eq. (Ill, 


at the lowest order in q the macroscopic electric field tends to a direction-dependent 
constant, 




q (q-Z*) 


K0 


eoD q' 


(26) 


where is defined in analogy with Eq. (131 and is the purely electronic rela¬ 


tive permittivity tensor. The main physical consequence of this is the well-known 
frequency splitting between longitudinal optical (TO) and transverse optical (TO) 


phonons in polar crystals. In particular, due to the contribution of Eq. (26) to the 


dynamical matrix, the LO dispersion curves behave nonanalytically already at zero 
order in q; that is, the eigenvalue and eigenvector associated with an LO branch gen¬ 
erally depends on the direction along which one approaches T. Such a nonanaliticity 
propagates directly to the electronic and lattice response functions described in the 
previous Section, and needs to be adequately treated in order to be able to apply the 
Taylor expansions described in Eq. (10) and (14).|^ 


® The response to an acoustic phonon in a nonpiezoelectric insulator is nonanalytic only at second 
order in q, so the situation appears here, at first sight, less serious than in the case of optical 
phonons. Recall, however, that the flexoelectric tensor is precisely an 0{q^) property, and therefore 
it is directly affected by such issues. 
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There are several ways to approach this problem. For example, the theory of 
Ref. [T^ was developed for purely transverse and longitudinal phonons separately, 
leading to flexoelectric coefficients defined at fixed E and D (electric displacement 
field) respectively. Here, we take the approach of removing the macroscopic E-fields 
in a physically meaningful way by assuming, following Martin,!^ that a very low 
density of free carriers is present in the insulating crystal, and that these are allowed 
to redistribute adiabatically in response to a phonon perturbation. In particular, 
within the Thomas-Fermi approximation, we write the free-carrier density as 

=-eofc^F^(r)> (27) 


where V{r) is as usual the electrostatic potential, and we suppose that the Fermi 
wavevector /ctf is much smaller than any reciprocal lattice vector of the crystal. In 
such a regime, the ground-state charge density and wavefunctions are essentially un¬ 
affected by the additional screening provided by the free-electron gas. Conversely, in 
the long-wave limit, the presence of the free carriers drastically alters the electrostat¬ 
ics; for example, the field of Eq. (26) becomeJi^ 


E 


q-i-o ^ q (q ■ Z*)k/3 

eoH -I- q • er • q 


(28) 


Such a modification has the following effects: 


• The macroscopic electric fields, and hence all the response properties of the 
crystal, become analytic functions of q. 

• The macroscopic electric fields vanish at zero and first order in q, and also 
at second order in q provided that we are considering an acoustic phonon 
branch. 

• Both the piezoelectric and flexoelectric tensors calculated in the presence 
of the free carrier gas are independent of fexFi and therefore can be un¬ 
ambiguously interpreted as the short-circuit versions of the corresponding 
electromechanical response functions. 


In the first-principles calculations, this is done in practice by simply suppressing the 
G = 0 contribution to the electrostatic energy when computing the self-consistent 
linear response; this has the same effect as introducing a low-density electron gas as 
described above. 

Based on the above discussion, it would be tempting to conclude that the flexoelec¬ 
tric tensor, like the piezoelectric tensor, is well defined under short-circuit electrical 
boundary conditions. In writing down Eq. (27), however, we assumed a particular 
type of carriers, namely electrons (not holes), and moreover that the band edge for 
those carriers, the conduction-band minimum (CBM), tracks with the macroscopic 
electrostatic potential of the crystal. In general, however, the CBM energy may shift 
relative to the local macroscopic potential as a result of a strain gradient, via the 
so-called deformation-potential effect. Thus, we can obtain a different flexoelectric 

^It is desirable to remove the macroscopic fields not only for practical reasons, i.e., to make the afore¬ 
mentioned Taylor expansions possible, but also because electromechanical tensors are traditionally 
defined in short-circuit electrical boundary conditions. 






Flexoelectricity 


9 


tensor depending on what band feature (CBM or other) we choose as the energy 
reference. We shall come back to this point in Sec. 

For a given energy reference, the bulk flexoelectric tensor fi is well defined in 
short-circuit (fixed E) boundary conditions. If fixed D boundary conditions are 
imposed along a specific direction q, the induced electric field (defined as the tilt of 
the corresponding reference potential) can be then easily calculated as|^ 




q 9a MaA,/37 

eo q • er • q 


(29) 


Note that Eq. (29) cannot be written in tensorial form, except for the simplest 
case of crystals with cubic symmetry, where the denominator reduces to a direction- 
independent constant. 


2.4. Lattice response 

To gain some insight into the nature of the lattice-mediated flexoelectric effect it 
is necessary to understand, in broad terms, the physics behind the internal-strain 
response (as described by the tensors N or L) to a strain gradient deformation. To 
that end, suppose that we perform a computational experiment where we statically 
freeze in a lattice distortion that corresponds to an acoustic]^ phonon truncated to 
first order in q, i.e., to the uniform-strain level, 

nL = (^a/3 + . (30) 

(In the simplest crystal structures, where the E tensor identically vanishes, this cor¬ 
responds to a purely elastic wave.) As we have perturbed the crystal from its equi¬ 
librium configuration, each atom in the lattice (identified, as usual, by a cell index 
I and a basis index k) will experience a restoring force f^.^. If the amplitude of the 
deformation is small (linear-response regime), such forces can be described, as usual, 
by a lattice-periodic (i.e., ^independent) function that is modulated by a complex 
phase with the same wavevector q as the perturbation. For small q, it can be shown 
that the magnitude of the induced forces scales as 0[q^) (first-order terms cannot be 
present, as we have assumed that uniform-strain effects are already included), and 
can be written as 


(31) 

Here is, by construction, the type-I flexoelectric force-response tensor. (The 


detailed derivation can be found in Ref. [I3ll4] .) 


Now one would be tempted, in close analogy to the piezoelectric case, to define the 
internal-strain response tensor N by means of the following linear system of equations. 


(T,(0) yy-ft' J_ rpK 
^ kock'^ q;/3,7A5 


(32) 


^Strictly speaking, this is the contribution from bulk effects; one cannot exclude surface contributions 

to the internal field, as we shall see in the later sections. _ 

^We assume that the long-range Coulomb fields have been removed; see Sec. 2.3 


for details. 











10 


M. Stengel and D. Vanderbilt 


where is the zone-center force-constant matrix.ji] Unfortunately, the above 

system is generally not solvable: the sublattice- (k-) sum of the T-tensor does not 
vanish, and the matrix is singular. (It is always characterized by three null 
eigenvalues, corresponding to rigid translations of the crystal as a whole.) As negative 
as it sounds, this is nonetheless an important result: it tells us that the internal-strain 
response to a static strain-gradient deformation is generally ill-defined. (We shall see 
later on that there are notable exceptions to this statement, though.) 

To understand what went wrong, let’s start all over again, but instead of consid¬ 
ering a static (frozen-in) deformation, take a dynamical one, i.e., a phonon mode. By 
performing a long-wave expansion of the equations of motion one obtainsl ^^ l ^'^ lfor the 
second-order eigendisplacements. 




P^'^P/3,7A 


rpK X ' rpK 

Z.^^a/3,7A) 

k' 


(33) 


where are atomic masses and M = all respect analogous 

to Eq. (321, except for the additional term that appears on the right-hand side (rhs) 
of the latter. It is trivial to check that the sublattice sum of the rhs now correctly 
vanishes, providing us with well-defined values (modulo a rigid translation) for the N- 
tensor components. This confirms our earlier suspicions that, unlike piezoelectricity, 
flexoelectricity is a genuinely dynamical effect: only in a sound wave are the internal 
strains well defined, and these internal strains depend explicitly on atomic masses. 
In retrospect, this conclusion is not entirely surprising. A uniform strain can always 
be generated and sustained by applying an appropriate distribution of external loads 
to the surface of the sample. This is not the case for a strain gradient: in general, 
a uniform force field applied to each material point of the sample is necessary to 
generate a given component of e/ 37 , a- Such a uniform force can be, e.g., generated 
by a gravitational fielcP^ or, as in the above example of the sound wave, by the 
acceleration of each material point during its periodic oscillation.l^^ In either case, 
the result directly depends on the atomic masses. 

To gain further insight into the physical nature of the mass-dependent term in 
Eq. (33), it is useful to write the same equation in type-II form. 


T K,' _ nr 

^KOiK,'p^pX^^-y ^aA,/37 


aX^^-y- 


(34) 


Here, C"* is the type-II flexoelectric force-response tensor, linked to T via the usual 
permutation of indices. 


C K _ rpK. I rpKs rpK 

ocX,0'y ^al3,yX ^a.'y,X^ -^OiX,j3'y 


(35) 


and Ca\,p-y is the macroscopic elastic tensor. To write Eq. (341 we have made use of 
the result 


E^«A./37 = ^C„y/3^, (36) 


i$(0) is the q ^ 0 limit of the matrix , , which is essentially a dynamical matrix with the mass 
prefactors set to unity. As for other quantities, is defined at vanishing macroscopic electric 

field, i.e., closed-circuit boundary conditions, appropriate for computing transverse optical phonon 
frequencies. 
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which directly relates flexoelectricity to elasticity^^ To justify such a sum rule recall 
that, in the context of linear elasticity, the stress tensor (which we allow to be 
inhomogeneous in space) is directly related to the elastic and strain tensors via 


^ —f^a/37A ^7A (l*) ■ 


(37) 


Recall also that the divergence of the stress tensor integrated over a finite region of 
space yields the net force acting on the corresponding volume element of the material. 


fa= d^rVpaap{r), 

Jn 


(38) 


By assuming that the crystal is homogeneous (i.e., that the elastic tensor is a con¬ 
stant), and by assuming that the deformation varies slowly over the volume of a 
primitive cell, we have 


, /ct(l') — ^CaP~f\£-f\,piT^ 


(39) 


Assuming that the force on individual atoms is exclusively produced by strain- 
gradient effects (which is justified, as the relaxations due to the local strain are 
already included), we can replace /” with the definition of the flexoelectric force- 


response tensor, and easily recover Eq. (36). Thus, in a hand-waving way, one can 


say that the type-II flexoelectric force-response tensor is a “sublattice-resolved” ver¬ 
sion of the macroscopic elastic coefficients. 

The dynamical nature of the flexoelectric tensor is worrisome if we are to use this 
theory to rationalize typical experiments - these are typically performed statically. 
As we shall see in the following, this is not a real issue. In a material is at static 
equilibrium there might be nonvanishing stress fields due to the application of external 
loads; nevertheless, the force acting on a material point must vanish everywhere in 
space. This leads to the following condition on the strain-gradient field, 


A,/97 ^/57,A (l*) — 6- 


(40) 


This means that two or more strain-gradient components will typically be present in 
any inhomogeneous strain field, in such a way that their respective net forces mutu¬ 


ally cancel. By using Eq. (401 it is straightforward to see that the mass dependence 


disappears from the resulting polarization field (as obtained by multiplying the flexo¬ 
electric tensor by the local strain-gradient tensor), confirming the internal consistency 
of the theory. 

The important message here is that, at the static level, we can define a number of 
effective flexoelectric coefficients; each of them will correspond to a linearly indepen¬ 


dent set of strain-gradient components that satishes Eq. (40). (An explicit example 


is provided in Sec. §) It is easy to see that the number of such effective static coefh- 
cients is always smaller than the number of independent components of the /.t-tensor. 
This means that the latter contains, in fact, more information than is actually needed 
to predict the outcome of a static measurement. This also means that, in order to 
determine the full flexoelectric tensor, one cannot rely on static experiments only; 
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additional dynamical data need to be combined with the static results!^ The result¬ 
ing values of the tensor components are always inherently dynamic quantities, even 
if static data are, in part, used to compute them. 


2.5. Electronic response 


While the lattice-mediated response has a straightforward physical interpretation 
(i.e., in terms of a polar distortion of the basis atoms that is induced by the macro¬ 
scopic strain gradient), the purely electronic response (given by the tensor /Sj) 
is far less intuitive, and therefore deserves a separate discussion. First, recall that 
fiaX Pj defined in terms of the second-order P-tensor, ■ To understand the 

physical meaning of the latter, consider the microscopic current density Ja{r) that 
is adiabatically induced when displacing an isolated atom {I, k) with velocity in 
the Cartesian direction ^)I2E3 


'Pa,K.p{Y) 


i9JQ,(r-f R;«,) 


diii 


K0 


(41) 


Provided that the macroscopic electric fields have been appropriately screened,^ one 
can introduc^i^ the moments of the vector field Va,Kp{j) at an arbitrary order n. 


^K7i.-7C ^ J d^rVa,Kp{r)rj, ... (42) 

Then, one can shovff^ that the resulting J-tensors coincide with the P-tensors of the 
same order apart from a trivial factor of volume, 


j: 


(n,7i...7„) ^ ^p(n.7l-7. 


a.,K,(3 


(43) 


This result tells us that the “frozen-ion” (in the sense specified in Ref. 12) contribu¬ 
tions to the piezoelectric and flexoelectric tensors are given in terms of the first and 
second moments of the current-density response to atomic displacements, respectively. 

Direct calculation of the P-tensors is technically challenging at the time of writing 
- the required current-density response functions are presently not available in the 
existing implementations of DFPT. To avoid this complication altogether, RestaP^^ 
proposed to determine the frozen-ion flexoelectric tensor via the sole knowledge of the 
charge-density response to an acoustic phonon, in close analogy with Martin’s classic 
treatment of the piezoelectric problem.l^ In particular, for an elemental crystal (this 
result was later generalized to arbitrary crystals by Hong and VanderbiltP^ Resta 
demonstrated that the longitudinal component of the response to a longitudinal strain 
gradient is given by 


p. = —Q^^\ 


(44) 


Here q indicates the spatial direction of interest, and indicates the correspond¬ 
ing third moment of the charge-density response to atomic displacement (dynamical 
octupole). 
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To derive this result in the context of the formalism of Sec. |2.2[ it is useful to 
introduce the charge-density response to the monochromatic lattice perturbation of 
Eq. 


—q • —(1,7) 


_(2,7A) 

2 


.q^qxqs _(3,7A5) 

■ * P.p 


(45) 


where the overline symbol implies cell averaging as in Eq. (13), and we have pushed 
the expansion up to third order in q. (The zero-order term vanishes because of the 
condition of charge conservation.) The p tensors are trivially related via 


(n,7i...7„) _ 1 ^(n,7i...7„) 

PkP — 


to the moments 


^(n,7i...7„) 

'^K.P 




= / d^rf^p{v)r. 


7i 


(46) 


(47) 


of the charge-density response function (r) ,|^ defined as the change in charge den¬ 
sity resulting from a single ionic displacement k/ 3. One can also show that the P- 
tensors and p-tensors are related bjff^ 


(n,7i...7jv) _ p("-b7l-..[7i]---7^) 

PkP A^-^7i,kP 

I 


(n > 1 ), 


(48) 


where the symbol [ 7 ;] indicates the absence of the element I in the list. Then one 
immediately has, for n = 3, 


72 . 7 A) .(2.a7) T(2.Aa) 

’^a,KP "'A,k/3 “I" '^j,kP 


_ ^(3,ajX) 
— ^kP 


(49) 


By applying Eq. (18) and Eq. (49) to the case of a longitudinal strain gradient oriented 


along q, one easily recovers Eq. (441. 


Unfortunately, it is not possible to invert Eq. ( |49[ ) and extract all components of 
the J*^^)-tensor from the octupolar response tensor, (The fact that contains 

more information than can be already appreciated by counting the maximum 
number of independent entries in either tensor: 54 in the former, 30 in the latter.Il^ 
Therefore, working only with the charge-density response at the bulk level is not a 
viable route to achieving full information over the frozen-ion (electronic) flexoelectric 
tensor, fi. 

Such a limitation can be circumvented, at least in cubic crystals, by considering 
a more general class of deformations that cannot be straightforwardly described as 
bulk acoustic phonons. For example, as we shall see in the next Section, the open- 
circuit internal field that is linearly induced by bending a free-standing slab is a 
bulk property of the material. Since the electric field is uniquely determined by the 
induced charge density this gives us, in principle, access to the transverse component 
of the electronic flexoelectric tensor without the need for calculating the polarization 
response. A bending deformation can be conveniently simulated (although at the 
price of a significantly higher computational cost) by adopting a slab geometry, and by 


J The function / is, in all respects, analogous to that introduced by Martin in his seminal work on 
piezoelectricity 0 
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performing a long-wave analysis analogous to that described here to the corresponding 
slab supercell. 

A formal derivation clarifying whether such a procedure does indeed yield the same 
transverse flexoelectric component as the P-response theory is still missing, due to 
subtleties at both the conceptual and technical levels. We shall refer to these issues 
again in the discussion following Eq. (101). In the remainder of the Chapter we shall 
disregard such issues, and provisionally assume that this relationship holds, i.e., that 
the bending-induced open-circuit (OC) E-field and the corresponding component of 
the flexoelectric tensor are related by 


t^xx,yy 


de,, 


(50) 


^yy,x 

[for a beam bent as in Fig. |^b)], where is the (isotropic) relative permittivity of 
the material. We shall use Eq. ( [50| from now on, whenever necessary, to resolve the 
aforementioned indeterminacy in the transverse components of /2. 

Note that this issue does not apply to the simpler case of the piezoelectric response. 
In fact, one can write that 


t(1^7) I t(1,«) 


= Q\ 


(51) 


(Recall that the basis sum of the tensors essentially coincides with the frozen-ion 
piezoelectric tensor, and that is the dynamical quadrupole tensor.) The above 
equation can be readily inverted. 


t(1.7) 

'^aP 



(2,q;7) 

(9 



(52) 


which provides an alternative derivation of Martin’s theorji^ of piezoelectricity. For 
completeness, it is useful to mention that, at order zero, the relationship between J- 
and Q-tensors is even more direct. 


t(0) _ ^(1,0:) _ 


(53) 


where Z* is the Born effective charge tensor associated with the k sublattice. Thus, 
both and can be regarded as higher-order generalizations of the dynamical 

charge concept, although starting from n = 2 (which is relevant for flexoelectricity) 
the former quantities generally carry more information than the latter ones. 


2.5.1. Spherical term, pseudopotential dependenee, and the noninteracting 
spherical-atom paradox 

As an illustration of the above derivations, it is useful in this context to work out 
a simple toy model that can be solved analytically; this will be also useful to point 
out some unconventional aspects of the flexoelectric response that have no connter- 
part in earlier theories of electromechanical effects in solids. We consider a rocksalt 
ionic crystal such as NaCl or MgO, and suppose that a longitudinal strain gradient 
develops along the (100) direction. Here we shall focus on electronic effects only, 
so that the atomic x coordinates undergo displacements that are a predetermined 
quadratic function of x with no further relaxations. For the time being we shall also 
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(a) 



(b) 







E 


Fig. 1. Simplified sketch of the planar-averaged electron potential energy, —eV{x) (black curves) 
for (a) an undistorted crystal, and (b) a crystal with a uniform longitudinal strain gradient. The 
red dashed lines show the macroscopic averages of the aforementioned functions. 


assume that the crystal is perfectly ionic, i.e., that its electronic charge density can 
be approximated by a superposition of spherical closed-shell ions whose shape is not 
altered by changes in bond distances, etc. With the above assumptions in mind, one 
can perform an average of the electrostatic potential in the yz planes, and express 
the result as a one-dimensional function of x. The atomic planes will appear as a 
periodic arrangement of potential wells (each well corresponding to a single charge- 
neutral monolayer), whose shape will reflect the radial distribution of electrons in the 
constituent ionic species. For the present purposes, the fine details of the potential 
wells are irrelevant; the only important quantity will be 

/ OO 

VyrL{x)dx, (54) 

-OO 

where is the j/z-averaged electric (Hartree) potential 14iL(r) generated by one 

monolayer, and —eVMh{x) is the corresponding electron potential energy. Thus, for 
purposes of illustration we can represent the potential-energy wells as nonoverlapping 
rectangular dips of area \K\, whose shape is fixed and independent of the surrounding 
neighbors, as sketched in Fig.[^a). As the wells are all identical and their separations 
are uniform in the undistorted crystal, the macroscopic electron potential energy 
obtained by convoluting the corresponding microscopic function with an appropriate 
low-pass filter,E^ shown as a dashed line in Fig. Ba). is constant. After freezing 
in the strain gradient deformation pattern, as shown in Fig. Bb), the interlayer 
distance increases linearly along the chosen axis, leading to a constant slope in the 
macroscopic electrostatic potential and, hence, to a uniform electric field throughout 
the bulk crystal. This result points to a nonzero flexoelectric coefhcient of purely 
electronic origin, since we explicitly neglected possible internal strains. 

There are several important questions that naturally arise at this stage. The first 
obvious one is whether (and if yes, how) the outcome of Fig. can be rationalized in 
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the context of the theory developed in this Chapter. A second and less obvious issue 
arises in pseudopotential-based first-principles calculations, where one may wonder 
whether and how the results depend on choice of pseudopotential. The third question 
concerns the physical nature of the electric field that we describe in Fig. [^b). Does 
it, for example, produce a direct force on charged particles such as electron and hole 
carriers and ionic cores? 

To answer the first question, it suffices to suppose that the electrostatic potential 
wells are generated by a regular lattice of spherical charge distributions. To make 
things simple, consider a monatomic lattice, as for a rare-gas solid, that we construct 
by periodically repeating a spherical charge distribution po(r). We assume that the 
volume of the unit cell D(r) depends smoothly on space as a result of an inhomoge¬ 


neous macroscopic deformation. One can show (see Supplementary Note 1 of Ref. 15) 


that the resulting macroscopically averaged (in three dimensions) electric potential 
is given by 

1 /■ ., o . . 1 


l/(r) = -- 


d'^s s'^pois) = 


rOr 


(55) 


6eon(r) J 6eoD(r) 

where Ol = 47 r / ds po{s) is the isotropic quadrupole moment|^of the static charge 
distribution po(r). Equivalently, it is the longitudinal component Ol = 


of the dynamical octupole tensor defined in Eq. (471, as follows from straightforward 
algebra. 

In the linear regime (small deformations) we have 


f2(r) ~ 0(1 -I- det [e(r)]). 


(56) 


which leads to the variation in the macroscopic electrostatic potential induced by the 
deformation, 


AE(r) = 


GenO 


det [e(r)]OL 


(57) 


Assuming that the crystal has cubic symmetry and making use of Eq. (501, this yields, 
after some algebra, two of the three independent components of the flexoelectric 
tensor 

Maa./3/3 ■ 


(58) 


(In the case of a biatomic ionic crystal one simply needs to replace Ol with the 
sublattice sum of the dynamical octupoles of the individual atoms.) By using the 
relationship between J-tensors and Q-tensors discussed earlier in this Section, it is 
not difficult to deduce that the third component, with a ^ j3, must be zero. 

Summarizing the above, the three independent components (longitudinal, transverse 
and shear) in a rigid-sphere crystal read as 


^xx,xx 


Ol 

en’ 


f^xx^yy 


Ol 

6D’ 


f^xy^xy 


= 0 . 


(59) 


This demonstrates that the effect illustrated in Eig.[2is indeed a natural consequence 
of the theory developed in this Chapter. 

'‘That is, the trace of the 3x3 second-moment tensor. 
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We now turn to the second question, concerning the use of pseudopotentials in 
first-principles calculations, as discussed in Ref.jl^ One aspect of the pseudopotential 
approximation is the replacement of the all-electron charge density by a pseudo 

charge density p^^(r) in the core region of the atom. Since these charge densities are 
essentially rigid and spherically symmetric, the above considerations apply to them. 
As a result, to compensate for the use of the pseudopotential, one should add a “rigid 
core correction” 

=i7Tjds s" [pt^{r) - pf (r)] (60) 


to the longitudinal dynamic octopole of each atom k to recover the all-electron result. 
This propagates into a change and to a change of 

P-xx,xa: and fl^^x,yy (but not p}^y^xy) by 

This rigid-core correction is not small, and is not independent of the details of 
pseudopotential construction. Therefore, two different calculations of the bulk flexo- 
electric response cannot be directly compared unless this correction has been applied 
in both cases. Nevertheless, as long as the same pseudopotential is consistently used in 
the calculation, predictions of physical, experimentally measurable quantities should 
not be affected by this correction. In particular, we shall see in Sec. (2.6) than 
makes an equal and opposite contribution to the surface contribution. Because of 
this cancellation, the total (bulk and surface) flexovoltage response [see Eq. (61)] can 
be computed without the need for including this correction. 

The third question, regarding the physical nature of the resulting electric field, 
requires taking a closer look at some earlier works on the theory of absolute deforma¬ 
tion potentials]^^\^^\ (These can be regarded as the foundation of the modern theory 
of flexoelectricity, even if they were aimed at addressing a slightly different physical 
problem.) In a nutshell, if we wish to draw a band diagram of a crystal subjected to 
a strain-gradient deformation, knowledge of the macroscopic electrostatic field is not 
sufficient. Indeed, the relative location of the valence-band maximum or conduction- 
band minimum with respect to the electrostatic reference is itself a function of the 
local strain (via the so-called band-structure term), which implies that each band 
will “see” a different electric field. This means that one band edge may be perfectly 
flat, and the corresponding carriers feel no force whatsoever, even while the other 
band edge and/or the mean electrostatic potential can be strongly tilted.]^ In fact, 
even a metal subjected to a strain gradient will generally have a nonzero internal 
macroscopic electric field arising from a gradient in the mean electrostatic potential, 
although no current will flow. Thus, one should be careful not to interpret the macro¬ 
scopic electric field produced by the flexoelectric effect in a longitudinal acoustic wave 
as a “real” physical field; it is just the tilt of some arbitrary reference energy that 
may have little to do with the phenomenon of interest in a given specific case. Just as 
for the notion of a “flexoelectric field,” care must be used when speaking of “short- 
circuit” and “open-circuit” electrical boundary conditions, as these are ambiguous in 


'Recall that we work in the framework of Eq. j50| , i.e., we extract the flexoelectric tensor components 
from the induced electrostatic potential, rather than from the polarization response. 

“The tilt of the mean electrostatic potential will also depend on choice of pseudopotentials when 
these are employed, but the tilt of the valence and conduction band edges will not. 
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Conduction 

band 



•/7 + 

Valence 

band 


Fig. 2. Sketch of valence bands (VB) and conduction bands (CB) for a slab with a strain gradient 
across its width. Dashed line indicates the macroscopic electron potential energy V{x) = —e4>[x). 
Hole carriers feel no force because the VB maximum is flat, while electron carriers feel a force to the 
right, both of which are contrary to naive expectations based on the electric field pointing to the 
right in the interior of the slab. 

the nonperiodic strain-gradient world. 

In light of the above arguments, it is legitimate to wonder whether the bulk flex- 
oelectric effect is experimentally measurable at all. In fact, there are good reasons to 
believe that the tilt of the mean electrostatic potential does not provide a realistic 
description of the response - at least no more realistic than other reference energies 
(e.g., the conduction band bottom, or the valence band top, or the Fermi level). First, 
as we have argued above, the present theory yields a finite open-circuit “flexoelectric 
field” even in a metal, which is physically inconsistent. Second, if we go back to the 
example of the noninteracting spherical atoms, there are apparent inconsistencies as 
well: Since we have assumed that each potential well is independent of its environ¬ 
ment, its motion cannot, in principle, be detected by an electrode that is placed at 
the far-away surface of the sample - and yet, the bulk flexoelectric coefficients do not 
vanish. We have, therefore, a sort of paradoxical situation, where the presence of a 
macroscopic electric field inside the material is indisputable, but at the same time 
there cannot be any open-circuit voltage, because of the hypothesis of rigid poten¬ 
tial wells (which excludes long-range effects). To resolve these paradoxes, and place 
the present theory in the right context regarding experimental measurements, it is 
necessary to account for surface effects. We shall see how to do this in the following 
Sec. EH 

2.6. Surface effects 

Knowing whether a given physical property is sensitive to the details of the sample 
surfaces is a matter of central importance in condensed matter theory. In the majority 
of cases (e.g., piezoelectricity), surfaces typically start to matter only at small length 
scales, where they are responsible for deviations in the measured property from the 
corresponding bulk value. There are situations, however, where such a sensitivity to 
the crystal termination persists up to the macroscopic scale; flexoelectricity belongs 
to this category. In the present Section we shall elaborate on this statement from 
a heuristic point of view, which is anyway sufficient to illustrate the most relevant 
physical ideas. A more formal discussion, based on a microscopic theory of the 
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response to deformations, will be presented in Sec. |2.7| and Sec. |2.8[ 

In order to calculate the flexoelectric response of a finite object such as a slab 
it is appropriate to consider, rather than the induced macroscopic polarization, the 
open-circuit voltage At^ produced by the deformation.]^ We shall only focus, in 
the following, on contributions that tend to a finite constant in the limit of infinite 
thickness t, and introduce the flexovoltage coefficient. 




1 dAV 
lim - - 


t—>-oo t x 


(61) 


Recall that e/ 37 ,a = dep-y/drx is the gradient of the symmetric strain tensor along 
the Cartesian direction r\, and x indicates the direction normal to the surface.]^ For 
simplicity, here we shall also restrict our analysis to strain gradients of the type eaa,x, 
i.e., a diagonal (either longitudinal or transverse) component of the symmetric strain 
tensor that is linearly growing across the slab thickness. (These are sufficient to 
describe the bending of a free-standing slab; a more general analysis, including the 
shear component, is deferred to Sec. |2.8| ) We shall write the flexovoltage coefficient 
as a sum of bulk and surface-specihc contributions. 


^XX^OLOC — 


bulk 

XX,aa 




surf 

XX,aa 5 


(62) 


whose explicit forms will be derived in the following paragraphs. 


2.6.1. Electronic surface response 


First let us consider only the purely electronic (frozen-ion) response. Strain gradients 
of the type eaa,x are governed by Eq. (50); in our present notation this implies that 
the (open-circuit) uniform electric field that builds up in the interior of the slab as 
a consequence of the deformation is uniquely given in terms of the bulk flexoelectric 
coefficient of the material and its macroscopic dielectric constant by 

^^aa X frozen—ion 


(63) 


Here eo and Cr are the vacuum and relative permittivities, respectively, while is the 
type-II flexoelectric tensor; as before, we use the bar symbol to distinguish frozen-ion 
quantities from fully relaxed ones. Since the electric held is minus the derivative of 
the potential, the bulk internal held contribution to the overall open-circuit voltage 
is then proportional to t, leading to a hnite contribution to the overall hexovoltage 
coefficient that we identify with 


-bulk 
r XX,aa 


^^slab 


de 


aa,x 


frozen—ion 


l^xx,aa 


(64) 


indicate here by AV" the total potential step that builds up, as a consequence of the mechanical 
deformation, between the two vacuum regions located at either side of the slab. 

°In spite of its notation, should not be thought as a tensor. First, the surface contribution 

depends on the specific details of the crystal termination, and is therefore not a simple function of the 
surface plane orientation. Second, the bulk contribution is defined in fixed-D boundary conditions 
and therefore it has a nonanalytic behavior [see Eq. ( |29[ |] in all materials except those characterized 
by cubic crystal symmetry. 
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Fig. 3. (a-b): Linear response of the electric field (a) and the electric potential (b) to a uniform 

strain applied to a film of thickness t. The function (where U indicates “uniform”) can be 
regarded as representing a surface piezoelectric response (surface dipole layer appearing in response 
to strain), (c-d): Linear response of the field (c) and the potential (d) in the case of a strain gradient. 
The sketch shows the derivative with respect to {eaa,xt) for a strain variation eaa{a:) = xeaa,x in 
a slab extending from —t/2 < x < t/2. The field response contains two contributions. The first is 
given by , appropriately scaled by the linearly varying local strain. Note that the induced surface 
dipoles, schematically illustrated by arrows, now point in the same direction. The other is given 
by genuine strain-gradient effects contained in E^, which essentially reflects the bulk fiexovoltage 
response. The resulting potential response in (d) thus consists in a macroscopic internal field plus a 
surface dipole contribution. 


The surface contribution in Eq. (62) originates from the fact that a surface 
can always be characterized by a potential offset (j) between the macroscopic potential 
just inside and just outside the surface, and that this offset is different for the two 
surfaces in the presence of a strain gradient. Consider first the case of a uniform 
strain Saa applied to a slab of thickness t, as shown in Fig. a-b). The figure shows 
the derivative of the macroscopic electric field (panel a) and electron potential energy 
(panel b) with respect to the applied uniform strain Saa , and is the corresponding 
derivative of the potential offset tp. The variation of (p with strain can be regarded as 
resulting from the fact that the surface, by virtue of its lack of inversion symmetry, 
is locally piezoelectric.]^ For the slab as a whole, however, a uniform strain does not 
produce a net voltage, since the induced potential offsets on either side of the slab 
cancel each other, consistent with the fact that the overall slab is nonpiezoelectric. 

In the case of a strain-gradient deformation, on the other hand, the local strains 
at the opposite surfaces are opposite in sign, and do not cancel out, as illustrated 


P In another language, we are basically describing a strain dependence of the surface work function, 
although technically the latter is referenced to the valence band maximum rather than the average 
potential in the subsurface region. 
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in Fig. [^c-d). The slab is taken to extend over —t/2 < x < t/2 with local strain 
£aa{x) = xeaa,x, reaching values of Saa = ^{t/2)eaa,x at the two surfaces. The 
figure shows the derivative of the field (panel c) and potential (panel d) with respect 
to {Saa,xt)- This means that the induced potential offsets at the two opposite surfaces 
have the same sign and add up in a flexoelectric experiment, leading to a surface 
contribution of the form 


-surf 
r xx,aoc 


d(j) 




frozen—ion 


The total flexovoltage coefficient then reads as 


^XX^OLOC 


l^xx,ota. 

eo^r 


—surf 

^XX^OLOL ' 


(65) 


( 66 ) 


The above derivation allows us to solve the paradoxes that we mentioned at the 
end of the previous Section. First, recall that we encountered some difficulties in 
giving a physical interpretation to the “internal electric field” that is induced by a 
strain gradient, as such a field depends on the reference energy (i.e., Bloch electrons 
in different eigenstates do not experience the same electrical force). This issue is 
easily solved by observing that the surface potential offset (j) suffers from the same 
ambiguity as the bulk flexoelectric field; we defined it relative to the macroscopically 
averaged electrostatic potential under the surface, but we could have used the valence 
or conduction band edge instead. It is easy to see that the respective ambiguities 
contained in the surface and bulk terms exactly cancel, yielding an overall flexovoltage 
coefficient that is uniquely defined. Next, we have observed that there is an apparent 
physical inconsistency in the rigid-spherical-atom model, in that there should be no 
overall voltage response to a strain gradient, and yet the bulk flexoelectric coefficient 
does not vanish. It is easy to see that, once the surface contribution is taken into 
account, the total flexovoltage response of a slab made of noninteracting spherical 
atom is zero as it should be. Indeed, when such a slab is subjected to a uniform 
strain, its surface potential voltage response is 


-surf -surf 

rxx^xx rxx,yy 


Ol 

beofl ’ 


(67) 


since a positive longitudinal or transverse strain increases the spacing between the 
atomic spheres and thereby reduces the surface potential offset. But, using Eqs. (58) 
and (64) (and the fact that Cr = 1 for this model), this is exactly leading 

to the claimed cancellation in Eq. (66). This cancellation also explains why the 
replacement of the all-electron by the pseudo core charge density in the context of 
pseudopotential calculations has no effect on the total flexovoltage response, so that 


the rigid-core correction of Eq. (60) can be neglected, as was claimed in Sec. 2.5.1 


Note that the spherical atom model, in spite of its simplicity, is crucial to un¬ 
derstand how flexoelectricity works in real materials. As we shall see in the results 
section, there generally tends to be a large cancellation between surface and bulk 
contributions to the flexoelectric effect. This happens because, even in covalently- 
bonded materials, the electronic charge distribution that is dragged along by each 
atom during its motion is largely constituted by a spherical shell, with comparatively 
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smaller aspherical components. Spherical objects do not contribute to the overall 
flexovoltage coefficient of a slab, hence the aforementioned cancellation. 

This gives a measure of the importance of the surface contribution - only when it 
is correctly taken into account together with the bulk term we obtain a meaningful 
physical quantity. Therefore, asking whether the surface contribution is “large or 
small” compared to the bulk effect is a poorly formulated question; the two must 
always go hand in hand. Instead, a more physically meaningful question is “How 
strong is the dependence of the surface contribution on its atomic and electronic 
structure?” 

Based on these considerations, one can attempt to give an answer to a long¬ 
standing question that has been somewhat controversial in recent years: “Is flex- 
oelectricity a bulk property?” As we said above, if by “flexoelectricity” we refer 
to the result of a typical flexoelectric experiment (i.e., where the induced current 
upon bending a short-circuited slab is measured), the answer is no. Conversely, if 
by the same name we call the current flowing through the bulk of the material while 
well-defined internal electrical boundary conditions are imposed, then the answer is 
yes. The problem is that the internal electrical boundary conditions depend on the 
externally-applied ones in a way that is surface-dependent, and unlike in the case of 
most known material properties, such a dependence persists in the limit of a macro- 
scopically thick sample. All in all, in the present context we would rather stay away 
from the traditional rigid classification into bulk properties and surface properties, 
as flexoelectricity, strictly speaking, does not belong cleanly to either category. 

2.6.2. Lattice surface response 

We now discuss how the above conclusions need to be modified when full ionic re¬ 
laxations are incorporated - these are, of course, of the utmost importance for a 
realistic description of the flexoelectric effect. Essentially, the above conclusion still 
hold, except for two important details: (i) the frozen-ion quantities (flexoelectric 
coefficient, dielectric constant, surface potential response) need to be replaced with 
their relaxed-ion counterparts; (ii) an effective deformation, given by an appropriate 
linear combination of, e.g., a longitudinal and transverse strain gradient, need to be 
considered in place of the individual tensor components. 

To illustrate the implications of (i) and (ii) in a practical situation, it is useful to 
work out the explicit formulas for the simplest case of an unsupported slab subjected 
to bending.]^ Linear elasticity dictates that a transverse strain gradient (correspond¬ 
ing to a “frozen-ion” bending deformation) at static equilibrium must be accompanied 
by a longitudinal strain gradient, which for most materials will have opposite sign 
compared to the transverse one. In fact, the top layers of the slab (“top” here means 
furthest from the bending center) are under tensile strain, and this typically induces 
a longitudinal contraction of such layers, whose magnitude is related to the Poisson’s 
ratio of the material. Conversely, the bottom layers are transversely compressed, 

^We shall exclusively focus, for the time being, on the p/ate-bending regime, where any deformation 
(e.g., anticlastic bending) along the main bending axis is forbidden. More general situations will be 
considered in the later Sections. 
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and will therefore expand longitudinally by an equal amount. This means that, to 
calculate the static flexovoltage coefficient of a bent slab, we need to consider the 
“effective” deformation 


^yy,x ^xx,x 

rather than the individual strain-gradient tensor components, where 

C 

^XX.VV 


c 

^x 


( 68 ) 


(69) 


is uniquely given by the elastic constants of the bulk material. Consequently, when 
the ions are relaxed, we shall be concerned with an effective flexovoltage coefficient 
reflecting the aforementioned mechanical equilibrium condition, 


^xx,ef[ — 


r^xx^efi 

eo^r 


dcj) 

de, 


eff 


where 




= 


(70) 


(71) 


refers to an analogous linear combination of the uniform strain components. 

The fact that, even at the level of the surface contribution, we have an effective 
response to a combined transverse and longitudinal strain is fully consistent with the 
behavior of an unsupported slab subjected to uniform in-plane tension. In such a 
situation, the relaxation will affect not only the surface atoms, but will also extend 
to the entire slab, leading to a contraction in the third dimension proportional to the 
bulk coefficient v. Thus, for a free film in a relaxed-ion context, it is only meaningful 
to consider the response of the surface potential offset (f to ees, and not to the 
individual £yy or e^x components; the former is precisely the quantity that enters the 


total flexovoltage coefficient in Eq. (701. 


Of course, one generally needs to consider more realistic mechanical boundary 
conditions than that of a free-standing film. In such cases, some of the specihcs of 
the above example are no longer valid (e.g., the absence of surface loads). Still, the 
points (i) and (ii) are applicable to the most general case. 


2.7. Electronic and lattice response revisited: Curvilinear coordinates 

In the early Sections of this Chapter we have described a fundamental theory of the 
bulk flexoelectric effect, based on a first-principles quantum-mechanical description 
of the insulating crystal. Later, in Sec. |2.6| we have argued, based on heuristic argu¬ 
ments, that there are important surface contributions to the flexoelectric response of 
a finite sample, and that these need to be accounted for when discussing experimental 
results. Here, we shall put the derivations of Sec. |2.6| on firmer theoretical grounds 
by developing an alternative approach. In particular, we shall clarify how to describe 
the microscopic charge and current responses to an arbitrary inhomogeneous strain 
field in terms of cell-periodic response functions. Such a formalism is necessary in 
order to treat, in full generality, the response of a finite (and hence, spatially inhomo¬ 
geneous) body to a deformation. As we shall see later, this will be useful not only for 
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Fig. 4. (a) Unperturbed charge density of model ID crystal composed of Gaussian charge packets, 

(b) Change of charge density in linear response to a uniform strain (solid line) or strain gradient 
(dashed line). 

the formal derivation of the surface contributions to the flexoelectric effect in finite 
samples, but also for the practical calculation of the transverse bulk components of 
the flexoelectric tensor. (Recall that such components are presently difficult to access 
at the bulk level.) Given its rather technical character, and the fact that the most 
relevant physical results have already been presented in Section |2.6[ this Section and 
the following can be skipped on a first reading. 

2.7.1. A simple one-dimensional example 

In order to establish a microscopic theory of deformations, the first issue one needs 
to address concerns the proper representation of the scalar and vector fields that 
describe the physical property of interest (e.g., atomic positions, electronic charge 
density, etc.). To appreciate the nature of the problem, it is useful to analyze the 
charge-density response of a simple lattice to a macroscopic deformation. Consider a 
one-dimensional chain of equally spaced atoms, which we represent as a regular array 
of Gaussian charge distributions as in Figurea). Its unperturbed charge density is 

p{^)='^Po{x- Rn), po{x) =, (72) 

where = na is an integer multiple of the lattice parameter a. Now we apply a 
uniform expansion to the chain by displacing each atom by 

Urfi eRji^ ('^^) 

and we look at how the charge density responds to such a perturbation . In the 
linear limit (small lambda) we obtain the response function dp{x)/d£ that is plotted 
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as the black curve in Fig. Bb)- The form of dp{x)/de is manifestly problematic: 
such a function grows linearly when moving away from the origin, i.e., it is clearly 
nonperiodic, which contrasts with the fact the system remains periodic after the 
application of the perturbation. Moreover, it introduces an undesirable dependence 
of the result on the arbitrary location of the coordinate origin. Such issues become 
even more severe when considering a strain-gradient perturbation of the type 

Un = \rI- (74) 

The charge density response, plotted as the red curve in Fig.|^b), now grows quadrat- 
ically with the value of the unperturbed atomic position, and extracting any relevant 
physical information from such a function appears difficult. 

The solution of the above problems comes from the realization that the fixed 
laboratory frame is a poor choice of coordinate system if we wish to represent the 
response to a macroscopic elastic deformation. In such a frame, the boundary atoms 
in a large crystallite have to move very far from their original location even if the 
applied strain is small; if we naively take the difference in the charge density from the 
original to the current state we obtain a result that has little physical meaning, and 
most likely will strongly deviate from the linear regime that we have in mind. A viable 
alternative is to treat an elastic deformation as a deformation of space, rather than 
an atomic displacement pattern. This implies applying a coordinate transformation 
that exactly reproduces the macroscopic elastic deformation.]^ From this viewpoint, 
the atoms do not explicitly move from their original location, although they do move 
with respect to the laboratory frame because the coordinate system itself is changing. 

To be explicit, consider a distortion r' = r -b u(r) that maps point r in the 
original periodic crystal into point r' of the distorted crystal, and such that a nucleus 
at R/k would be carried to = Rj^ + u(RiK) if one neglects the additional internal 
displacements arising from the lattice effects described in Sec. |2.4| If the initial charge 
density poijc) were also carried along by this distortion, the new charge density would 
be 

Prefir') = Pair) det"\h) (75) 

where the Jacobian factor involving = dr'^/drp = 5ap -\- dua/drp is needed to 
reflect the dilution or concentration of charge density. In fact, the actual charge 
density pir') has to be computed from the appropriate physical laws (e.g., first- 
principles DFT calculations), so it will not be equal to pref(r^)- However, we may 
hope that the difference p(r') — pi.ef(rO small, and we want to express this difference 
in terms of the original spatial variable r. This is conveniently done by defining 

p{r) = p(r') det(h) (76) 

so that our small quantity is Ap(r) = pir) — po{r). Note that /3(r) describes the actual 
charge density after the deformation, but transformed back to the original coordinate 
system; the hat symbol is used henceforth to highlight quantities that describe the 
transformed system from the curvilinear-coordinate point of view. 


"^Recall that a deformation of a continuum is a 3D-3D mapping of each material point to its perturbed 
location, i.e., it has the exact same mathematical form as a coordinate transformation. 
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Fig. 5. (a) Same as in Fig. |^a). (b) Change in charge density, when expressed in transformed 

coordinates, for a uniform strain (solid black line) or a strain gradient (dashed dark gray line). The 
latter, while not periodic, can be expressed as a sum of black and light gray contributions (the latter 
was magnified by a factor of 50), as explained in the text. 

In Fig. we again perform the same analysis as in Fig. illustrating how the 
use of coordinate transformations effectively solves the problems that we pointed out 
earlier. Panel (a) shows the same charge density at rest. As before, in this model 
we assume that the actual charge densities shift rigidly with the nuclei. In panel 
(b) we plot as the black curve the induced density dp{x)/ds for a uniform strain. 
The response is now periodic and much smaller in magnitude than before (note the 
scale change). We shall denote this response function as p^{x), where ‘U’ indicates a 
‘uniform’ strain. The response to a strain gradient, shown as the dashed red curve, is 
still not periodic, although it now has a milder dependence on the spatial coordinate, 
growing only linearly rather than quadratically with x. Remarkably, however, we can 
write this response as a linear combination of two cell-periodic functions, 

^^^^=xp^{x)-^p^{x), (77) 

where p^(x) is the same as above (response to uniform strain), and p^{x) is a new 
quantity, reflecting the genuine strain-gradient effects (shown as a thick light gray 
curve in the figure, where it has been magnified by a factor of 50 to better illustrate 
its functional form). Since we are considering a uniform strain gradient above, we 
have e{x) = xp, so that we can write 

Ap{x) = e{x)p^{x) + + • ■ • (78) 

In other words, we have achieved a closed expression for the induced charge density 
Ap(x) that depends only on proper measures of the local deformation, with the only 
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hypothesis that the local strain £{x) varies slowly on the scale of the interatomic 
spacings. 

Several questions naturally arise from the above discussion. First, how general is 
such an analysis? For our illustrative example above we have used a trivially simple 
system, and a single (longitudinal) strain (or strain gradient) type, so it is legitimate 
to wonder whether the same procedure is applicable to a full first-principles simulation 
in 3D. Second, what do we do with Ap{x) once we have calculated it? To make the 
discussion relevant for flexoelectricity it is necessary to trace a direct link between 
Ap{x) and measurable electrical quantities, such as the macroscopic polarization in 
short circuit, or the induced voltage in open circuit. We shall address both questions 
in the remainder of this Section. 


2.7.2. General formalism in three dimensions 

Regarding the general applicability of the coordinate transformation method, there 
are several conceivable ways to proceed. One could, for example, directly incorporate 
the curvilinear-coordinates formalism at the level of the Kohn-Sham equations (bor¬ 
rowing from the adaptive coordinate scheme of GygP^ and, in a similar spirit as in 
Ref. |21[ directly perform the perturbation expansion with respect to the metric tensor 
and its gradients. Alternatively - and we shall follow this latter strategy throughout 
this Chapter - one can go back to the phonon analysis that we have introduced in 
Sec. |2.2[ this time focusing on the microscopic charge-density response functions; the 
challenge here lies in converting these to the curvilinear representation outlined in 
this Section. We thus consider a deformation 


r'pir) =rp + Upe^' 

which generates a simple frozen phonon 


.p = Upe 




(79) 


(80) 


For the moment we neglect the internal displacements leading to the lattice response 
of Sec. 2.4 so that Eq. (80) is equivalent to Eqs. (9p0 ) with the q-dependent terms 
neglected, but they will be restored shortly in Sec. 2.7.4| 

In the linear limit, the charge density responds as 


p(r) = po{r) + p^ir) e^^ ’^ 


(81) 


where the cell-periodic part Pp{r) gets modulated by the same phase factor as in 
Eq. (79). Inserting this in Eq. (IT^ gives 


(82) 


/5(r) = (po(r') -H C/; 3 p^(r')e*‘i '’') (l -t 


where the last term in parentheses is the value of det(h) resulting from Eq. ( [7^ . We 
now expand the cell-periodic response function up to second order in q. 




1 . 7 )/ 


Q-yQX (2,7A) 

-^Pp 


2 


(83) 










28 


M. Stengel and D. Vanderbilt 


Since we are only collecting terms to first order in U in Eq. (82), we can ignore the 
distinction between r and r' in the cross terms, but for the direct term we have 


^0(1-') = ^0(1- + Ue*'! '’) = po(r) - C// 3 p^°Hr)e 


(84) 


where we have used that dppo{r) = —p^^\r). Collecting all the terms linear in U, we 
obtain 


A/3(r) 


iqppiy) - iq^p^p''^\ 


r) 


979a (2,7A) 

-^Pp 



(85) 


The p^^^ terms have now canceled, as expected from the fact that the coordinate 
transformation has removed the translational part from the response. 

After observing that the unsymmetrized strain and strain gradient are related to 
partial derivatives of the displacement held. 


£p^{r) = iq^Upe^^'", ( 86 ) 

??/3.7A(r) = -979AC//3e*'‘■^ (87) 


we can readily write 

A/5(r) = £p^{r) \6p^p{r) - p^^’^^(r) 


(2.7A) . X 

n Pp 


( 88 ) 


Finally, one can replace ip^ with the symmetrized counterpart, £p^ (the quantity 
in the square brackets is invariant upon /Sy exchang^i^, and replace 9 / 3 , 7 A with the 
type-II strain gradient tensor e/ 37 ,a- This leads to an expression that is in all respects 
analogous to Eq. (78), 

Ap(r) = £p^{v)p^p^[r) + 

where the uniform (U) and gradient (G) terms are dehned as follows. 


PM^t^) = ^PiP{^) - P^p'^'\^) 




(2.A/3) 

7 


(r) -pf’^^^(r) 


(89) 

(90) 

(91) 


This result formalizes and generalizes the arguments of the hrst part of this Section: it 
shows that the microscopic charge density response to an arbitrary inhomogeneous de¬ 
formation can indeed be computed (and rigorously expressed in terms of well-dehned 
response quantities) in a Hrst-principles context, and for an arbitrary combination of 
the relevant 3D deformation tensor components. 

As a side note, one can show that the quantity Pp^{v) essentially coincides (apart 
from a trivial scaling factor) with the Hrst-order charge-density as defined by Hamann, 
Wu, Rabe and VanderbiH!^ within their linear-response theory of strain based on the 
metric tensor. It will be interesting in the near future to draw even closer connections 
between the two formalisms, which bear several similarities at the conceptual level. 
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2.7.3. Microscopic polarization response 


In the above derivations we have focused on the charge-density response of the sys¬ 
tem to an inhomogeneous deformation, but we could have worked just as well with 
the microscopic polarization response instead. This quantity is well-defined only for 
infinitesimal transformations, otherwise it depends on the specific path followed by 
the system during its evolution; this is not an issue here, since we are exclusively 
interested in the linear-response regime. The microscopic polarization P(r) is re¬ 
lated to the adiabatic current-density response of the system to a time-dependent 
perturbation. Thus, in order to construct an appropriate definition of this quantity 
in a generic curvilinear frame, we need first to examine the transformation laws of 
the current-density field J(r). To this end, let r' = r -|- u(r,f) be a generic time- 
dependent coordinate transformation, which we suppose to coincide, as usual, with 
the displacement field associated with the mechanical deformation of the sample. (We 
suppose now that such a deformation happens slowly over a finite interval of time.) 
In a curvilinear framework the four-current, defined as = {p, ji, j 2 , js), transforms 
as 


J'" = 


dx'^ 

dx°‘ 


J“ defi 


dx^ 

dx'y 


(92) 


where x^ = (t,xi,X 2 ,x^) is the coordinate four-vector and the barred (unbarred) 
symbols refer to the deformed (original) frame. We work here in the nonrelativistic 
limit with t = t and r independent of t, so that 


p = pdet ^(h), (93) 

Ji = (p^ +hijJj^ det“^(h). (94) 


(Latin indices refer to the three-dimensional Cartesian space.) This leads to the 
definitions 


p = p det(h), 

i, = (h-i)., 


- _duj 


det(h). 


(95) 

(96) 


The above expression for the curvilinear-frame charge density p coincides with that 


postulated earlier in Eq. (76), showing that this definition is, in fact, dictated by 


the fundamental transformation laws of a scalar density field. Equation (96), on the 


other hand, gives the desired expression for the curvilinear-frame current density Ji, 
which we will use in the following to derive the microscopic polarization response to 
an acoustic phonon perturbation. 

As in the former case of the charge density response, we consider a monochromatic 


acoustic phonon as in Eq. (79), again without internal cell relaxations. This time, 


however, we allow the amplitude of the displacement to depend on time. 


r' = r + u{r,t), up{r,t) = Up{t)e"^''‘ 


( 97 ) 
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(Henceforth we shall go back to using Greek indices for three-dimensional space co¬ 
ordinates, as we did in the previous sections.) In the adiabatic limit, one has 

9P(r) 


J(r,t) = Ui3{t)- 


dUp ’ 


and, after dropping all terms that are quadratic in either Up or Up [this implies 
setting hij = 5ij in Eq. (96)], 

J(r,t) = Up{t)^^^-tJ{t)poir). 

Now, the microscopic polarization response can be written, as usual, as a cell-periodic 
part times a phase, P(r) = Up P]g(r), which immediately leads to 


Pair) = UpU^ '^ -Pa,/ 3(0 - ^aPPoir) 


(98) 


By following the same steps as for the case of the charge-density response, we now 
proceed to expand P^ pir) in powers of q, 

P^(r) = P^°)(r) - *g^P^'’^^(r) - 


>ir). 


(99) 


From translational invariance, it is then easy to show that the zeroth-order term 

Pa,pir) = Sappoir), (100) 


exactly cancels with the last term involving po in Eq. (98). Eventually, we arrive at 
a provisional result for the linearly induced polarization currents in the curvilinear 
frame of the deformed body in the form 


APa = ep^ir)P^^p^{r) + 
where the cell-periodic vector fields P^P are 


P 


aX,0'y 


(r), 


Pa,p^ir) = -Piy\r), 


P 


aA,/37 


(r) = ^ + PiV^Hr) - P^T^^) 


2 1 


( 101 ) 

( 102 ) 

(103) 


To arrive at this equation, however, we have had to assume that the currents generated 
by a global rotation are the same as those obtained by rigidly rotating a classical 
charge density that is eq ual to the true quantum-mechanical one. This assumption, 
which is implicit in Eq. (50), was used to conclude that PiV\r) = Pa’p\r), and 
hence to replace the unsymmetrized (e) with the symmetrized (e) strain tensor (see 
Sec. V.C of Ref.[^. While such an assumption was indeed valid in the charge-density 
case, it is not obvious that it is justifiable in the present case of the microscopic 
polarization current. Discussing this point in detail would take us far from the scope 
of the present Chapter; nevertheless, the reader is warned that there are still some 
unresolved formal issues in the theory of the current-density response. 

Regardless of such issues, one can show that the functions Ap and AP enjoy the 
fundamental relationship 


V • AP(r) = —A/5(r), 


(104) 
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where V • A indicates the divergence of the vector field A in the curvilinear frame. 
(The hat is used to emphasize that the differentiation is with respect to r rather than 
r'.) This reflects a well-known fact, which is important in the specific context of flex¬ 
oelectricity: the induced charge density can be readily deduced from the polarization, 
but not the other way around. As a matter of fact, taking the divergence annihilates 
the solenoidal part of the AP-field, which does contribute to the bulk flexoelectric 
tensor. (Recall the relationship between dynamical octupoles and second moments of 
the current-density response: the additional information contained in the latter can 
indeed be ascribed to divergenceless polarization currents that arise in response to 
an atomic displacement.) 


2.7.4. Atomic relaxations 


We now return to the inclusion of the internal atomic relaxations, and describe how 
they can be conveniently incorporated in the aforementioned theory; the practical 
implications regarding surface effects will be discussed in Sec. |2.8[ 

The first important observation is that, given a strain field ep^{v,t) that depends 
slowly on space and time, the internal-strain tensors that were introduced in Sec¬ 
tion |2.2| can readily be identified with the microscopic lattice response of the crystal 
in the curvilinear frame of the deformed body, 

4 + ■ ■ ■ ( 105 ) 

In particular, for a macroscopic strain gradient, the above relationship reduces to 


dul, 


de 


Pi A 


— Rikx -i- L 


-'a.\,0'y ’ 


(106) 


Next, it is easy to show that Eq. (101) still holds, provided that we replace the purely 
electronic response functions with their relaxed-ion counterparts. 


( 107 ) 

pXpi^^) = P^x.Pii^) - Pi].''Ji^Kpi + Pi°U^)PpX,Pi- (108) 

(We have used the bar symbol here to denote the purely electronic response 

functions and thereby distinguish them from the fully relaxed quantities.) This for¬ 
mally extends the microscopic linear-response theory discussed in this Section to the 
relaxed-ion case. 


2.7.5. Electrostatics in a curved space 

Having established a convenient form for the microscopic charge and polarization 
response functions, we still need to figure out how to use them, e.g., how to calculate 
the voltage response Viv). This requires some attention, given the fact that we are 
no longer working in the Cartesian frame of the laboratory. In particular, one needs 
to replace Gauss’s law with its “curvilinear” generalizatiorESESl 

V • (e • E) = p. 


(109) 
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where the vacuum permittivity has been replaced with the tensor 

e = eodet(h)g"\ (110) 

the hat on V is again a reminder that the gradient is in the curvilinear frame, g = 
h • h"* *" is the metric of the deformation, and 

^a(r) = KpEp{v') (111) 

is the transformed electric field.|^ As Vy(r) = —E(r), the above strategy allows 
one to compute the induced electrostatic potential from the induced polarization (or 
charge density). 

In the linear limit of small deformations, one has 

eoV • (AE + AE“®‘) = A;5(r) (112) 

where AE(r) = — V[Ay(r)] is minus the (curvilinear) gradient of the induced elec¬ 
trostatic potential, AI^, and AE™®*, coming from the linearization of e, reads as 

AA“®‘(r) = e/37(r) [5p^Ea{v) - SapE^iv) - 5^aEp{v )], (113) 


where E^{v) is the electric field in the unperturbed system. The choice of notation 
AE“e* is meant to suggest a “metric contribution to the curvilinear-frame electric 
field,” but this is somewhat problematic as it is not an irrotational field. Alternatively, 
cqAE™®* could be regarded as a “metric contribution” to the polarization, since one 
can rewrite Eq. (112) in terms of AP as 

eoV • AE =-V • (AP-f eoAE™®*) . (114) 


However, this is not entirely satisfactory either, as AE™®* does not really originate 
from the displacement of charged particles. It is probably most appropriate to inter¬ 
pret cqAE™®* as a displacement current arising from the effective change of permit¬ 
tivity associated with the deformation of the reference frame. 

In any case, Eq. (112) shows that the induced potential AE(r) contains, in addi¬ 
tion to contributions from the rearrangement of the electron cloud occurring during 
the deformation (these are contained in Ap), also a term that depends on the local 
variation of the metric at fixed charge density. We shall come back to this point in the 
discussion of surface contributions to the flexoelectric effect in SectionNote that, 
by construction, the microscopic electric field response to an arbitrary deformation 
enjoys an analogous representation as the charge density response, 

dep^(v) 




(r) = ep^{r)E'^^p^{v) ■ 


drx 


-E 


aA,/37 


(r), 


(115) 


where both and are lattice-periodic functions whose explicit expressions can 
be readily derived from Eq. (112). 


® One can arrive at Eq. | |109| l by observing that, in a curvilinear frame, Poisson’s equation reads as 
y/ 9 ~^ (^y/ 99 ^^di,V^ = —pjeQ where g = det (g) = det^ (h), i.e., the Laplacian must be replaced 

with the Laplace-Beltrami operator. Then, by defining p = -^pp, Su = —d^V^ and = €. 0 y/gg^^, 
one immediately recovers Eq. ( |109[ |. 

*It is interesting to note the close connection between the curvilinear-frame electrical qu antit ies (E 
and P) described here and the reduced electrical variables (respectively Sa and pa) of Ref. |24[ where 
a linear mapping was implicitly used to connect lattice and Cartesian coordinates. 
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2.7.6. Treatment of the macroscopic electric fields 


Eq. (112) specifies AE(r) modulo an r-independent integration constant, AE whose 
value is fixed by the electrical boundary conditions (EEC) of the problem. The 
electronic response functions are typically defined (and calculated) by assuming AE = 
0, i.e., short-circuit (SC) EBCs, but any other EEC choice can be recovered if the 
charge-density (and/or polarization) response to a macroscopic electric field is known. 

For example, in the case of the polarization one can write 


AP„(r) = AP“^(r) + AEpP^^ir), 


(116) 


where Pa'^ (r) = dPa (r)/9P/3 is the microscopic P response to an applied field along 
yjlHl contribution of the macroscopic field can be readily incorporated into the 
strain-gradient term (P® in this case), as the macroscopic electric field response in 
a nonpiezoelectric crystal vanishes at the uniform-strain level. Therefore Eqs. (89), 
(101), and (115) remain valid in arbitrary EEC. 


2.8. Surface effects in curvilinear coordinates 

Recall that, in order to quantify the flexoelectric response of a free-standing slab, 
we have introduced the flexovoltage coefficient <Px\,i 3 -y Here we shall demonstrate 
how this quantity can be rigorously derived in the context of the microscopic theory 
developed in the previous Section. To express ^xX,M terms of well-defined response 
functions of the system, we shall follow the strategy of Sec. |2.7[ now specializing to 
the case of a supercell geometry. As before, we shall derive the microscopic response 
of the system (charge density, polarization, and atomic displacements) to a strain- 
gradient deformation via a long-wave analysis of its acoustic phonons. With the help 
of a coordinate transformation to the curvilinear frame of the perturbed body, one 
can express such microscopic response functions in terms of “proper” measures of the 
local deformation, i.e., in a translationally and rotationally invariant form, 

A/(r) = ef,^{r)f^^{r) -t + ■ ■ ■ (117) 

Here / can stand for the charge density (p), polarization (P) or electric field (E) 
expressed in the curvilinear frame. Note that the cell-periodic functions f^ and 
/®, referring respectively to uniform and gradient terms, are characterized by an 
oscillatory behavior on the scale of an interatomic distance, due to the discreteness 
of the atomic lattice. As it is customary in the space-resolved analysis of many 
other physical properties (e.g., dielectric response), we shall assume in the following 
(unless otherwise specified) that such oscillations have been filtered out by means of 
an appropriate nanosmoothin^HEH technique. 
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Fig. 6. Schematic representation of the three types of strain gradient described in the text, (a) 
Longitudinal, exx,x- (b) Transverse, eyy^x- (c) Shear, £xy,y The x and y axes correspond to the 
horizontal and vertical directions in the figure respectively. (Adapted from Ref. |27[) 


2.8.1. Surface polarization and metric 


To illustrate the above arguments in the present context, consider a symmetrically 
terminated slab of a cubic material (we assume that the surfaces are parallel to the yz 
plane), and perturb it with a strain-gradient deformation. We assume for the moment 
that the ionic coordinates simply follow the deformation; we shall lift this limitation 
in the next subsection. In order to calculate the flexovoltage coefficient of the slab, we 
shall first derive the electric field E(r) induced by the deformation under open-circuit 
electrical boundary conditions. Then, by performing a line integral of E(r) across 
the slab thickness, one can readily obtain the desired value of <Pqa,/ 37 - To calculate 
E(r) we need, in turn, two basic ingredients: the microscopic polarization response, 
AP(r), and the “metric” contribution to the polarization, AE™®*. Regarding the 
former, after nanosmoothing are functions of x only, and we can write 


dPo^jr) 


frozen—ion 




(118) 


The two response functions and Pa\,pj{^) have the physical meaning of a 

local piezoelectric and flexoelectric coefficient, respectively. Note that Pa,pj{x) differs 
from zero only in the vicinity of the surface, as the bulk material is nonpiezoelectric. 
The metric term, on the other hand, reads as 




CX {dpyEa(^tc') Sappy {x'j SyaPpi^x)^ ■ 


(119) 


This quantity, just like P^, is active only at the surface: the electric field of the 
undistorted slab is nonzero (and directed perpendicular to the surface plane) only in 
a small region where the crystal lattice is perturbed by the truncation of the bonding 
network. 

The details of the derivation differ, from now on, depending on the specific type 
of flexovoltage response that one wishes to calculate. It can be shown that, given 

^Note that AE is first-order in the perturbation, and therefore its contribution to the electronic 
response functions is only due to the microscopic dielectric properties of the unperturbed system. 
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Fig. 7. Schematic representation of the polarization fields induced by different macroscopic de¬ 
formations of a slab of thickness t, drawn in the undistorted reference frame. The x coordinate is 
vertical here. Left panels (a-b) illustrate uniform strains; right panels (c-d) refer to uniform strain 
gradients. Top (a,c) are transverse deformations (the situation is qualitatively identical in the lon¬ 
gitudinal case, not shown), bottom (b,d) are shear patterns. The surface region and corresponding 
polarization field are indicated by light gray shading and black arrows. White arrows on a dark gray 
background refer to the bulk region. The dashed black frames indicate the type of deformation in 
each case. (Adapted from Ref. |15[ ) 


the symmetry of the slab, there are only three types of strain-gradient deformation 
that yield a net open-circuit voltage in a cubic material: a variation of €xx with x 
(longitudinal), Syy with x (transverse), or e^y with y (shear), which are responsible for 
the flexovoltages ^xx,xx^ ^xx.yy, ^^id y^xy.xy respectively, in the notation of Eq. (611. 

Longitudinal and transverse cases. In these two cases (which we shall indicate 
as XX, aa), the system remains periodic in-plane, and the problem becomes essentially 
one-dimensional. We have, in particular. 


dPx{x) 


FI 


De 

L/Cn 


dE^^\x) 


= xP, 


.u 

x,ctoc \ 


i^)+p: 

= xEx{x){l — 26, 


G 

XX,Q.OL 




( 120 ) 

( 121 ) 


(where FI is shorthand for ‘frozen-ion.’) The polarization response functions and 
xP^ P’^ are schematically illustrated for the transverse case in Figs. 0:^) and (c) 
respectively. Given that both vector fields are irrotational and vanish at infinity, one 


can safely simplify Eq. (114) by removing the divergence sign on both sides to get 
dExix) 


dSn 




1x1 -b 

Y V*^/ I XX,OLOi 


(x)] - xEx{x){l - 2Sax)- (122) 


The frozen-ion flexovoltage coefficient of Eq. (61) can then be calculated by writing 
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the open-circuit potential associated with the above field as 

/*+oo 

7G 


AV = - 




(123) 


where 


E^.^4x) = --P^,ac.ix) - E,{x){l - 2S^,), 
eo 

1 -G 


U 

x,aa 

G 


(x) = _ ix') 

xx,acx\'^J xx,aa\‘^J' 

eo 

The above functions enjoy a number of useful properties: 


(124) 

(125) 


(i) E^^^^ (x) vanishes everywhere except for a small region near the surface at 

X ~ Etl2\ 

(ii) E^^^^{x - t/2) = —X + t/2) is antisymmetric, and independent of t 

for a sufficiently thick slab; 

(iii) E2x^aaix) corresponds to minus the bulk fiexocoupling coefficient in the slab 
interior, 


Exx,acxi.X ^ 0) (p. 


bulk 
XX,act 


.bulk 


eoEr 


and only deviates from this value in a small region near the surface; 

(iv) E2x,aaix — t/2) = E2x,aai~x +1/2) is Symmetric, and again independent of 
t for a sufficiently thick slab. 

Based on these observations, in the limit of large slab thickness one can approximate 


Eq. (123) as 


AE - -t 


P + OO 

/ dx E 2 ^^{x) + tiplf^a- 

Jo 


(126) 


(We assume that a; = 0 is the center of the slab and x = -l-oo is deep in the vacuum 
region.) As the integral in the last equation is independent of t, we can readily write 


^+oo 




dx E2„t,ix) + 


(127) 


whence we obtain 


-surf 
r xx,occx 


^+oo 


dxE^ctaix). 


(128) 


The last equation states that the surface contribution to the fiexovoltage response 
of a slab corresponds to minus the line integral of the induced electric field upon 
application of a uniform strain. The latter is, of course, the electrostatic potential 


offset response to uniform strain, d/i/dsaa, that we already discussed in Sec. 2.6 


We have, therefore, rigorously demonstrated that the fiexovoltage response of the 
slab indeed contains both bulk and surface contributions, and that their nature is 


correctly described by Eq. (66). Note that the derivations presented here, in addition 
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to corroborating the arguments of Sec. |2.6[ allow us to make one step further and 
split iiito 9, polarization current and a metric term, 


d(j) 

L/Cq-ck 


1 

eo 







(129) 


In the last term on the right, 4>o = — dx Ex{x) is the potential offset before the 
perturbation. 

In summary, in the present case (longitudinal or transverse strain gradient) the 
induced electric field in the interior of the film is a bulk property of the material 
- it is given by the flexoelectric coefficient divided by the macroscopic dielectric 
constant. The surface contribution, on the other hand, acts as an induced potential 
offset that grows linearly with slab thickness, and therefore scales similarly to the 
bulk contribution, as illustrated in Fig. 

Shear case. The case of a shear deformation {xy, xy) is qualitatively different 
from the former two cases .Q Here we have 


dPJv) 


FI 


3e 

^^xy,y 

aF;“®*(r) 

3e 


— dayyPyxyiP) + Pc 


ay,xy 


{x), 


— dayyEx{p^^- 


(130) 

(131) 


Figures Hb) and (d) illustrate the polarization fields that are linearly induced by 
shear deformations. By taking the divergence of the above vector fields, we can write 
the curvilinear Poisson’s equation, Eq. (114), as a function of x only,|3 




cq- 


Ax) 


dx 


= -P,) 


dpG 


(x) 


dx 


+ ^Opx{x). 


(132) 


[We have used the short-hand notation E^y^xyix) = dEx{x)/dexy,y] Assuming that 
the field vanishes inside the slab,|2we can then calculate the macroscopic electric field 
in the vacuum region. 


eo- 


dEx{x = -boo) 


de 


xy,y 


FI, sc 


= cr. 


surf 


,bulk 


(133) 


(SC stands for ‘short-circuit.’) The three quantities on the right-hand side have the 
physical dimension of a surface charge density and are given by 


bulk _ -;bulk 
^xy,xy Fxy^xy^ 


surf 
^ xy,xy 



dxPy^xyix), 


met 

^xy,xy 


^+oo 


eo 


dx Ex(x) 


—eo<)>o- 


(134) 

(135) 

(136) 


"^We mention this case for completeness, as it is not relevant for an unsupported slab after full atomic 
relaxation, provided that we consider a region that lies far (compared to the thickness, t) from the 
edges and/or the mechanical loading points. We shall come back to this issue in Sec. 2.8.2 
'^Note that the partial derivatives along y of both and vanish identically, as these nanos- 
moothed functions are periodic in plane, and therefore only depend on x. 

^This is the natural choice for the electrical boundary conditions when considering shear strain 
gradients - recall that they directly relate to transverse acoustic phonons. 
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Fig. 8. Flexoelectric response of a slab of noninteracting spheres. The total induced potential 
is decomposed into three contributions, consistent with the formalism developed in the main text. 
In all panels, the vertical axis corresponds to the potential (dashed line indicates the zero), and 
the horizontal axis is the spatial coordinate x along the surface normal. The combined effect of 
the surface (including both induced polarization-charge and metric contributions) and the bulk 
flexoelectric response yields a vanishing bias potential, regardless of the type of strain gradient. 
(From Ref. |15| ) 


In order to derive the flexocoupling coefficient, we need to switch to open-circuit 
boundary conditions by imposing an external electric field that exactly cancels the 
above vacuum field. We obtain 


_ 1 

^xy,xy — Z“ 
eoCr 


f^xx.aa. 



dxPy^^yix) - eocpo 


(137) 


The total surface contribution coming from both the polarization currents and 
from the metric is thus 


-surf 

rxy,xy 


1 

Coer 



dxPy^^yix) - eo^!)o 


(138) 


Note that, in contrast with the transverse and longitudinal cases, the internal electric 
field is no longer a bulk property here; the surface terms contained in manifest 

themselves as surface charge densities that tend to a constant in the limit of large 
slab thickness t, rather than dipole densities that grow linearly with t. (In either 
case, the surface contribution to the flexovoltage response of the slab scales similarly 
to the bulk term for increasing t.) These, unlike in the previous two cases, need to 
be divided by the bulk permittivity, as the bulk material dielectrically screens the 
additional electric field produced by surface effects. 

Spherical atom model. The formalism that we have developed in this Section 
allows us to complete the solution of the toy model that we described at the end of 
Section 2.6 consisting of a finite slab made of a lattice of rigid (and noninteract¬ 
ing) closed-shell atoms. The solutions for all the contributions to the flexovoltage 
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response, now including the shear case and the aforementioned separation of the sur¬ 
face term into polarization charge and metric terms, are schematically illustrated in 
Fig. i. (The details of the derivations can be found in the Supplementary Notes 
of Ref. [^) As expected, in all cases the net voltage vanishes, consistent with the 
physical expectations (a rigid displacement of spherical charge distributions cannot 
lead to a long-range electrical perturbation). By the same token, the pseudopotential 
rigid-core correction of Eq. (601 has no effect on the net flexovoltage. Note, how¬ 


ever, that the bulk, surface-metric and surface-polarization terms cancel each other 
in a nontrivial way depending on the strain-gradient component, indicating that a 
consistent treatment of all three terms is crucially important for having a physically 
meaningful solution. 


2.8.2. Atomic relaxations 

The contribution of atomic relaxations to the flexovoltage coefficient of a bent slab has 


been extensively treated in Sec. 2.6.2 It is easy to show that, by using the formalism 
presented in Sec. 2.7.4 we recover Eq. (701, which describes the total response in 


terms of bulk- and surface-specific quantities. What remains to be discussed is the 
shear case. In Sec. |2.8.f] we postulated that this type of strain-gradient deformation 
is not relevant for a fully relaxed unsupported slab. Here we shall substantiate this 


statement in light of the results presented so far. Recall Eq. (106), which describes 


the microscopic atomic relaxation pattern induced by a strain gradient in terms of 
the internal-strain response tensors F and L, and let Xi^ and Yzk denote the x and y 
components of In the case of a shear strain gradient of the type £xy,y in Fig.|^c), 


Eq. (106) reads as 


de 


= Yu 


^y,v 


■pK. I T K 

^ ctxy ' -^ocy,xy' 


(139) 


Now, regardless of the microscopic details of the slab, rotational invariance dictates 
that 


Kxy = -XoJ, 


OK^ay ? 


(140) 


i.e., under a uniform shear the slab rigidly rotates to accommodate the deformation 
of the supercell, without feeling any restoring force because the repeated images of 


the slab are decoupled. By combining Eqs. (139) and (140) we obtain 

3£xv.v 


= + L" 


(141) 


Now, recall that a macroscopic strain gradient can be written in terms of the 
components of the type-I strain-gradient tensor as 




Eq. (141) states that a shear strain gradient of amplitude rix,yy = y is always accom¬ 


panied, in a fully relaxed unsupported film, by a second strain gradient component 
of the type rjy^xy = The overall effect, in type-II notation, is that of a negative 
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(a) __ (b) 



Fig. 9. Sketch of slab subjected to a periodic transverse strain of the type shown in Fig.[^b), or 
a negative shear strain of the type shown in Fig. j^c). After internal atomic relaxations, the two 
configurations become equivalent. 


transverse strain gradient, Syy^x = This means that, for a free-standing film, the 
shear case reduces exactly to the transverse one. The basic concept is illustrated in 
Fig.[9l where we conipare the configurations obtained by periodically subjecting a 
slab to a transverse strain gradient Syy^x as in Fig. [^b), or a (negative) shear strain 
gradient e^y^y of the kind shown in Fig. [^c). If internal atomic relaxations are al¬ 
lowed while still preserving the overall undulation along y, the two configurations will 
clearly relax to the exact same geometry. 

Thus, we have rigorously demonstrated the result that we heuristically presented 
in Sec. |2.6| fiexoelectric effects in a free-standing film of sufficiently high symme¬ 
try (e.g., cubic or in-plane hexagonal) are governed by only one response coefficient 
^xx,yy, as given by Eq. (701. The induced voltage at a given location is then given by 
^xx,yy{t/^y +t/^z), where ^y = = ^~zz,x are the radii of curvature (along 

the Cartesian axes) of the film at that specific point.|^ This includes the plate-bending 
and beam-bending limits as special cases. 

It is interesting to note that, in contrast with what happens in the bulk, here we 
have a notable case where the fiexoelectric effects induced by a sound wave are identi¬ 
cal to those associated with a static deformation. (Equivalently, one can say that the 
same strain gradient field can be induced either by dynamic or static means.) Indeed, 
any two-dimensional object such as a slab is characterized by a transverse acoustic 
phonon branch, usually referred to as ZA, with zero sound velocity, corresponding to 
a bending mode. A long-wavelength ZA phonon coincides, therefore, with the static 
bending case described above, and produces the same fiexoelectric response. 


2.9. Summary 

In this Section we have presented a fundamental theory of fiexoelectricity, based on 
a quantum-mechanical description of the electronic and lattice response to a strain- 
gradient perturbation. In particular, we have used a long-wave expansion of acoustic 

^One could equivalently choose different orthogonal axes, e.g., those corresponding to the princi¬ 
pal curvatures of the surface. Since 1/^y + is the trace of the shape operator, the result is 

independent of such a choice. 
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Fig. 10. Supercell models of the SrO- (a) and Ti02-terminated (b) SrTiOs slabs. Ti and O atoms 
are represented as white squares and gray circles respectively; Sr atoms are not shown. (Adapted 
from Ref. |27| ) 

phonons to derive, in the linear limit, the relevant electromechanical response func¬ 
tions of a crystalline solid. Our formalism is fully general, and correctly recovers 
earlier theories of piezoelectricity as a special case. 

In order to address some conceptual issues (e.g., regarding the role of the surfaces, 
or regarding the calculation of some components of the bulk flexoelectric tensor that 
are presently difficult to access) we have gone a step further, and developed a fully 
microscopic theory of the linear response to an inhomogeneous strain field. In this 
context, we have demonstrated that the use of curvilinear coordinate frames greatly 
facilitates the representation of the relevant physical fields and their response to 
mechanical deformation. The latter methodological tools are applicable well beyond 
the specifics of flexoelectricity, and may find application in related research areas, 
such as flexomagnetism.l^ 


3. Application to SrTiOa 


In this Section we shall demonstrate the theory developed so far by applying it to 
SrTiOa, one of the most important materials in the context of flexoelectricity, and 
the best known experimentally. In order to quantify the importance of surface effects, 
we shall consider a slab geometry, and two different lattice terminations (either of 


the SrO or Ti02 type), as illustrated in Fig. 10 


3.1. General methodology 

Our goal is to calculate the total flexovoltage response of either SrTiOa slab to a 
bending deformation in the limit of large thickness. We shall do this by taking into 
account the effect of full atomic relaxation, under the initial hypothesis that our slab 
behaves as a (We shall see in Sec. 3.3.3 that the beam-bending limit can 

easily be recovered by rescaling the plate-bending coefficient by a constant.) This 


®'This means that along the direction parallel to the bending axis the system is clamped (i.e., no 
anticlastic bending is allowed). 
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requires the combination of three different computational frameworks, as detailed in 
the following. 


3.1.1. Bulk calculations 


Here we perform a number of calculations on a primitive unit cell of bulk SrTiOa. 
This is primarily aimed at calculating the bulk flexoelectric tensor via a long-wave 
expansion of acoustic phonons. Acoustic phonons are treated at the linear-response 
level by means of density-functional perturbation theory as implemented in a modified 
version of the ABINIT packag^^ in which the contribution of the macroscopic electric 
fields has been removed according to the discussion of Sec. 2.3 In particular, We 


choose a small star of wavevectors q surrounding the T point in the Brillouin zone, 

q= —^(±1,0; ±1,0; 0), 

and perform a full linear-response calculation for each of these points. (In practice, we 
make full use of symmetries to minimize the number of actual calculations.) Next, 
we perform a long-wave expansion of the charge-density response and interatomic 
force constants and extract the second-order-in-q coefficients via numerical differ¬ 
entiation with respect to q. We obtain: (i) the flexoelectric force-response tensor 


p-i (311, from which and then the internal-strain tensor L'^x p-f 

are constructed via Eqs. ( 34|[35 l; and (ii) the charge-density response tensors 

and corresponding respectively to the Born effective charge tensor 

and the dynamical octupole tensor, via Eq. (45). 


By combining the internal-strain tensor with the Born effective charges one can 
readily obtain the lattice-mediated contributions to the flexoelectric tensor as ex¬ 
plained in Sec. |2.2| The octupole tensor, on the other hand, provides us with only 
partial information on the electronic (frozen-ion) flexoelectric tensor. In particular, 
only the longitudinal component of the electronic flexoelectric tensor, /Iq, along an 
arbitrary direction q can be inferred from the two linearly independent entries of 
■ Following Hong and Vanderbilt, (SI we define 

MLI = /I(IOO)) ilL2 = 2/2(110) ~ M(ioo)- 

These are related to the components of the type-H flexoelectric tensor, /rL 2 bjEH 


t^XX.XX MLIi (142) 

t^xx^yy ± ‘^l^xy,xy ML2- (143) 

Thus, in order to determine the transverse and shear components ft^xx,yy P^xy,xy 
independently, an additional calculation is necessary; this will be addressed shortly 
in Sec. 13.1.21 

In addition to the above calculations, which are based on the methodology de¬ 
scribed in this Chapter, we also need a bulk-level calculation of some auxiliary quan¬ 
tities by means of more established techniques. Specifically, we extract the high- 
frequency dielectric constant from a separate linear-response treatment of the 
electric-field perturbation. At the same time we obtain a redundant set of Z* p^ 
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tensor elements, which are useful for assessing the quality of the numerical differen¬ 
tiation at first order in q performed in Sec. |3.1.1| above. Similarly, we carry out an 
independent calculation of the elastic tensor Ca\,py via finite differences with respect 
to applied strain; this allows us to check the second-order-in-q calculations of the 
force-response tensors since these quantities are directly related by the sum 

rule in Eq. (ISbl). 


3.1.2. Truncated-bulk slab calculations 


Here we carry out calculations similar to those of Sec. 3.1.1 
supercell. 


but now on a slab 


This step is aimed at determining the transverse and shear components 
of the bulk electronic (frozen-ion) flexoelectric tensorIn fact, the two independent 
components of the bulk dynamical octupole tensor that we calculated above 

are not sufficient to determine the three independent entries of the bulk tensor. 

We are able to circumvent this limitation by recourse to a series of calculations on 
a slab geometry in which we determine the charge-density response, both in the 
bulk and at the surface, to longitudinal, transverse and shear strain gradients. A 
calculation of the flexoelectrically induced open-circuit electric field in the interior of 
the film, which relates [based on Eq. ( |50[ )] directly to the corresponding component 
of the bulk flexoelectric tensor in two cases out of three (longitudinal and transverse), 
allows us to obtain the missing component]^ of The key point here is that the 
missing divergence-free component of the induced polarization current, which is not 
currently available from bulk-level calculations, manifests itself as a surface charge 
density, whose influence is readily apparent in the slab supercell geometry. Note that 
the specifics of the surface structure should not matter in these calculations. Thus, 
we choose the geometry that ensures the best convergence of the inner open-circuit 
field as a function of slab thickness, i.e., a truncated-bulk structure. (We perform 
such an analysis on both SrO- and Ti02-terminated slabs, in order to verify that the 
results are indeed surface-independent as we expect.) 

In practice, we use the same star of q-points surrounding T as in the bulk calcu¬ 
lations described above. This time, however, we neglect the information on the force 
constants and only focus on the charge-density response of the system. We need to 
analyze such a response at the microscopic level, by using the curvilinear-coordinate 

Of the two relevant response functions, (x) 


formalism of Sec. 


2.7 


and Sec. 


2.8 


and p‘^{x), only the latter is really an issue, as p^{x) can be straightforwardly cal¬ 
culated as the response to a uniform strain.]^ The result of the second-order Taylor 
expansion in q yields p^{x), and this (together with p^) is then used to calculate the 


^ It may seem odd to use a slab supercell to calculate a bulk-specific quantity; this is indeed a 
temporary work-around, which will no longer be necessary once a proper theory of the current- 
density response becomes available. 

^Strictly speaking, only the transverse component is really needed, as the longitudinal component 
calculated in this way is redundant with the value that we already calculated at the bulk level. 
We shall use this as a test to assess the numerical accuracy of our calculations. 

^We calculated p^{x) separately by using standard ground-state calculations where we took finite 
differences in the strain. We found that this latter procedure yields slightly better accuracy than 
the long-wave method described above. 
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electric-field response functions 

Note, however, that due to the removal of the macroscopic electric fields in the 
phonon calculation£3IIIlZI (as required to perform the aforementioned Taylor expan¬ 
sions in q, see Sec. 2.3), short-circuit electrical boundary conditions are enforced by 
construction on the calculated and E'^. This means that there are nonvanishing 
macroscopic electric fields in both the vacuum and the slab interior, and these fields 
show an undesirable dependence on the supercell geometry (vacuum and slab thick¬ 
nesses). To have a physically well-defined (and geometry-independent) value of the 
internal field we need to enforce open-circuit electrical boundary conditions. We do 
this by applying an external field to the system that is exactly opposite to the cal¬ 
culated vacuum field. To determine the charge redistribution induced in the system 
upon application of an external field, we perform a separate linear-response calcula¬ 
tion of the local electric field response to a macroscopic electric displacement field D. 
This is nothing but the local inverse dielectric permittivity of the slab supercell, 

dEx{x) 1 

We then use 


E^’^^ix) = E^’^^ix) - Ef’^^i+^)-e-\x), 

where x = -foo corresponds, as usual, to the vacuum region. When referring to E^{x) 
in the following, we shall implicitly assume that we are speaking of the open-circuit 
version E^’^^ix). 


3.1.3. Relaxed-ion slab calculations 


Now that we have all the necessary bulk-specific information in hand, we still need 
to determine the surface-specific contributions to the flexovoltage coefficient 
We shall compute as the induced electrostatic potential offset upon application 
of a uniform effective strain {Syy = Seg] e^x = ~^£es) to a free-standing slab with 
(001) surface orientation. This quantity can be conveniently accessed by means of 
a standard plane-wave code; no linear-response features are needed. In particular, 
we take a slab supercell corresponding to a periodic lattice of alternating SrTiOa 
and vacuum layers, and first calculate the electronic and structural ground state by 
setting the in-plane lattice parameter to the equilibrium bulk value. We then apply 
a small positive or negative strain of the type 


e 



where e^g is a small dimensionless number, typically eeff(±) = ±0.001. (We find it 
computationally advantageous to preserve the fourfold axis of the SrTiOa surface by 
applying an isotropic in-plane strain.) 

^Recall that we need to consider, for a bent slab at mechanical equilibrium, an effective combination 
of transverse and longitudinal strain-gradient deformations, £yy,x = ^xx,x = —where 

V = Cyy^xxjC-xx,XX‘ 
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In each perturbed configuration, we first calculate the electronic ground state 
with the reduced coordinates of the atoms kept fixed to their unperturbed values; 
the resulting electrostatic potential profile is then processed by means of macroscopic 
averagin^HEU to extract the perturbed frozen-ion (FI) surface potential offsets. Next, 
we let the atoms relax to their new equilibrium positions in the strained lattices, and 
repeat the macroscopic averaging procedure to obtain the relaxed-ion (RI) offsets. 
Finally, we numerically differentiate the perturbed offsets (both FI and RI) to obtain 
their corresponding first-order variation, 

surf ^ - 4 >{-) 

^ ' 2|eeff| 

where <(>(±) refers to the surface potential offset at positive or negative strain.]^ This 
procedure readily yields the RI and FI values of The lattice-mediated (LM) 

values are simply calculated as the difference of the RI and FI ones. Of course, the 
slab needs to be sufficiently thick in order for the inner layers to be truly bulk-like, 
i.e., unaffected by the atomic distortions that originate from the surface truncation 
of the bonding network. 


3.2. Computational parameters 


We use the local-density approximatioiP^ to density-functional theory. The interac¬ 
tions between valence electrons and ionic cores are described by separable norm- 
conserving pseudopotentials in the Troullier-Martin^^ form, generated with the 
fhi98PP code.^ The and shells of Sr and Ti, respectively, are ex¬ 

plicitly treated as valence electrons. The reference states (numbers in parentheses 
indicate the core radius in bohr) of the isolated neutral atom used in the pseudopo¬ 
tential generation are 2s(1.4), 2p(1.4) and 3d(1.4) for O, 4s(1.5), 4p(1.5) and 4d(2.0) 
for Sr and 3s(1.3), 3p(1.3) and 3d(1.3) for Ti. The local angular-momentum channel 
is 1 = 2 for Sr and O and Z = 0 for Ti. The rigid-core corrections of Eq. (60) are not 
included in the presented results. The cutoff for the wavefunction plane-wave basis 
is set to 150 Ry in the slab calculations. (A test calculation with a 300 Ry cutoff did 
not show appreciable changes in the calculated electronic response functions; the 300 
Ry cutoff was, nonetheless, necessary to ensure satisfactory accuracy in the force- 
response tensor at the bulk level.) The surface Brillouin zone of the slab supercell 
is sampled by means of a 8 x 8 Monkhorst-Pack gridj^ for the bulk primitive cell 
we use a sampling of up to 12 x 12 x 12 fc-points. The finite-difference parameter in 
the long-wave expansion, q, is set to 0.01 (tests with q = 0.02 or q = 0.03 indicated a 
convergence better than 1% in the calculated electronic response functions; smaller 
values of q were found to yield less accurate results because of the excessive numer¬ 
ical noise). The lattice parameter of the cubic cell is set to ao=7.268 bohr, which 
corresponds to the calculated equilibrium value. 


^ As a technical note, many first-principles codes use the Ewald procedure to calculate the self- 
consistent electrostatic potential. This involves adding to the electronic density a lattice of spherical 
Gaussian compensating charges, whose spurious con trib ution must be removed from the calculated 
value of See the Supplementary Notes of Ref. |27|for details. 
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Table 1. Force-response tensor 
of bulk SrTiOa in short-circuit boundary 
conditions. Ol, 02 and 03 refer to oxygen 
atoms forming x-, y- or 2 :-oriented Ti-O-Ti 
bonds, respectively. All values are in eV. 


Atom 

{xx, xx) 

(xy,xy) 

{xx,yy) 

Sr 

-24.9 

7.9 

-28.7 

Ti 

-67.9 

3.8 

-102.3 

Ol 

159.3 

15.3 

97.4 

02 

35.2 

17.3 

42.3 

03 

35.2 

-0.9 

30.9 


Table 2. Summary of the linear-response data obtained from the 
long-wave (LW) approach at the bulk level, compared to the results 
of Hong and VanderbillBJ (HV) for the same quantities. Open-circuit 
electrical boundary conditions are enforced on the longitudinal response 
functions (LI and L2). The force response to a shear strain gradi¬ 
ent (S) is quoted in short circuit. The oxyge n m odes ^3 = xqi and 
^4 = (xq 2 + XQ;^)/y/2 are defined following Ref. |l4| is in V; other 

values are reported in eV. 



Ll(LW) 

Ll(HV) 

L2(LW) 

L2(HV) 

S(LW) 

S(HV) 

^bulk 

-16.15 

-16.25 

-18.07 

-18.17 

- 

- 

Sr 

16.3 

17.0 

33.2 

35.7 

7.9 

8.4 

Ti 

49.1 

52.3 

36.3 

38.9 

3.8 

3.0 

6 

67.2 

68.7 

24.9 

13.1 

15.3 

15.7 

^4 

3.0 

3.6 

22.7 

18.2 

11.6 

12.0 


The supercell models are based on the schematic illustrations of Fig. [Tola-b). For 
the truncated-bulk linear-response calculations we use 5.5-unit-cell (uc) thick SrTiOa 
slabs alternating with vacuum layers whose thickness is set to 2.5 uc. Of course, 
both (slab and vacuum) thicknesses are intended as convergence parameters in our 
calculations, whose scope is to describe the thermodynamic limit of a macroscopic 
slab. Tests with thinner slabs and thicker vacuum layers (up to 3.5 uc) showed optimal 
convergence for the aforementioned values of these parameters (again, better than 
1%). For the relaxed-ion slabs, we use 7.5uc-thick slabs with 3.5uc-thick vacuum 
layers. 


3.3. Results 

3.3.1. Bulk calculations 

In Table [2 we report the relevant values of the force-response tensor of bulk SrTiOa, 
calculated by using the long-wave method described in Sec. |2.2| In Tablej^we compare 
the above physical quantities to the analogous ones that were calculated in Ref.[^ To 
perform the comparison we first recast the force-response components into a tensorial 
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Table 3. Calculated Born effective charges and dielectric properties of bulk 

SrTiOs. 



7* 

Zxi 

^Ol 

7* 

^02 

•^03 

Er (static) 

2.5548 

7.2455 

-5.7027 

-2.0488 

-2.0488 6.1785 

1657 


Table 4. Calculated elastic tensor of bulk 
SrTiOs- The two rows refer to the bulk 
force-response calculation (“Force”) and to 
a direct bulk calculation where we took fi¬ 
nite differences of the calculated stress tensor 
while varying the strain around the equilib¬ 
rium cubic configuration (“Strain”). Values 
are in GPa. 


Method 

{xx, xx) 

{xy, xy) 

(xx, yy) 

Force 

385.3 

122.2 

111.7 

Strain 

386.2 

122.4 

112.6 


representation that follows the same prescriptions as Eq. (142) and (143), 


^L1 ^xx,xx ’ 

^L2 ^xx^yy ' xy,xy' 


(144) 

(145) 


Then, we convert the longitudinal quantities LI and L2 from fixed-E or short-circuit 
(SC) to fixed-D or open-circuit (OC) boundary conditions by using [see Eq. (106) of 


Ref.[^ 


Q(OC)=C£(SC)-t^t"“^Z:, (146) 

where Z* is the Born effective charge (calculated values are reported in Table , 
j^buik jg purely electronic flexovoltage coefficient, and L stands for either LI or 
L2. The calculated values of reported in Table for direct comparison 

to those reported by Hong and Vanderbilt.^ The agreement is overall very good, 
especially considering the different computational strategy, first-principles code and 
pseudopotentials that were used in Ref. 

As a numerical test of the calculated force-response tensor (Table[^, in Table|4] we 
report the elastic constants of bulk SrTiOa that we computed in two different ways: 
either as a first derivative of the stress with respect to the applied strain (“strain”) or 
by using the sum rule of Eq. ( [3^ (“force”). The agreement is excellent (better than 
1%), confirming the high numerical quality of the calculation. Note that the choice of 
the electrical boundary conditions is irrelevant for this test, as the sublattice sum of 
the atomic forces induced by a hypothetical electric field vanishes due to the acoustic 
sum rule. 
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Fig. 11. Electric field response to mechanical deformations. The (s.-b) and (c-d) response 
functions are shown for a SrO- (a,c) and Ti02-terminated (b,d) slab. Solid black, dashed black and 
solid gray curves refer to longitudinal, transverse and shear deformations, respectively. The location 
of the SrO (dashed) and Ti02 (solid) atomic layers is indicated by vertical lines (only half of the 
symmetric slab is shown). (Adapted from Ref. |27[ ) 


Table 5. Frozen-ion flexovoltage coefficients of a truncat¬ 
ed-bulk SrTiOa slab. To compute we used ¥^li^L 2 

as reported in Table [SI and E^^^yy = 15.08 V [extracted 
from Fig. [file -d)]. (Lj7(T) and (S) stands for longitudinal, 
transverse and shear, respectively. Units of Volts are used 
throughout. 



J^bulk 

V 

SrO 

surf 

Ti02 

(total) 

SrO Ti02 

XX, XX (L) 

-16.15 

14.36 

16.95 

-1.80 

0.80 

^^,yy (T) 

-15.08 

15.68 

12.45 

0.61 

-2.63 

xy,^y (S) 

-1.50 

-2.38 

-0.51 

-3.88 

-2.01 


3.3.2. Truncated-bulk slab calculations 


In Fig. ll'a-d) we plot the calculated corresponding to either a SrO- or 


a Ti02-terminated slab and to each of the three types of imposed strain gradients 
shown in Fig. (with no internal relaxations allowed). As anticipated in Sec. 2.8.1 


there is an important qualitative difference between the case of the longitudinal or 
transverse response, where the strain gradient is oriented along the surface normal, 
and that of the shear response, where it is directed in plane. 

In the former case, pp{x) (describing the E-field response to a uniform strain) 
yields the surface contribution to the flexovoltage coefficient, via Eq. (1281, 

1^ while the functions provide us with the sought-after information on the 


s Note, however, that here we are dealing with a truncated-bulk slab, whose surface atomic coor¬ 
dinates were artificially frozen to ideal bulk positions. The surface contributions that one extracts 
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bulk flexovoltage coefficient of SrTiOa, 

= ~Exx,0p{^ = 0 )- 


Note that the x) functions are roughly uniform inside the film, which in¬ 

dicates that the slab is thick enough to display bulk properties therein, and zero 
outside, consistent with the open-circuit electrical boundary conditions that were 
enforced. Moreover, the uniform internal held appears to be nicely independent of 
the surface termination for the longitudinal and transverse deformations, which is a 
further important consistency test for our computational approach. 

In the shear case, however, the hexoelectric held depends on both bulk and surface- 
specihc properties,^ and such a termination dependence is clear from a comparison 
of the gray curves in Fig. c) and (d). From the electric-held response functions 
of Fig. m c-d) we can thus only extract the total hexovoltage coefficient of the slab. 


(x = 0). To separate ipxy,xy into bulk and surface terms it suffices. 


^xy,xy — E^y^^y 

however, to complement the above data with the values that we calculated at 

the bulk level. Indeed, by replacing the hexoelectric tensor components in Eq. (142) 


and Eq. (1431 with the corresponding hexovoltage coefficients, we have 


, hulk _ ,>ulk 
t’IjI t’xx.xx'i 

, ^bulk , ^bulk I 

V^L2 ^xx,yy ^ 


hulk 


2(p 
r xy,xy 


(147) 

(148) 


Eq. (147) constitutes a useful consistency check of the methodology, as is redun- 

Eq. (1481, on the other hand, yields 
from the slab calculations. 


dant with the already calculated value of (p. 


,bulk 

XX,XX ' 


bunr 


^bulk 


the desired value of since we already know Pxx,yy 

Finally, we use Pxy,xy = -Efy^’xy to infer p'^xy^xy = ^xy,xy - ^xy,xy 

Our results for the bulk, surface, and total hexovoltage coefficients of the 
truncated-bulk, frozen-ion deformation of a SrTiOa slab are summarized in Table 
At the bulk level, it is interesting to note the relatively small magnitude of the shear 
coefficients and p’^xyxy compared to both the longitudinal and the transverse 

ones. Meanwhile, in the latter two cases there is a substantial cancellation between 
bulk and surface terms; as a result, the values of the total hexovoltage coefficients p 
are all comparable in magnitude. This fact can be rationalized by observing that the 
linear response to atomic displacements in a ionic (or partially ionic) solid is largely 
dominated by the rigid displacement of an approximately spherical charge density 
distribution surrounding each atom. The spherical contribution, which is typically 
large and negative,ti^ shows up in and with opposite sign in in the 

shear case neither the bulk nor the surface term are affected (see Sec. 2.5.1). Re¬ 
markably, the resulting values of p depend strongly on the details of the surface, and 
in some cases even have opposite signs in the SrO- and Ti02-terminated slabs. Such 
a conclusion, in fact, persists after we take into account the full relaxation of the 
atomic structure; we shall demonstrate this point in the following paragraphs. 
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Table 7. Flexovoltage coefficients of a relaxed 
SrTiOs slab. The frozen-ion (FI), lattice-mediated 
(LM) and total relaxed-ion (RI=FI-|-LM) values of 
the bulk, surface and total slab response are reported. 
Units of Volts are used throughout. 



^bulk 


surf 

tsp (total) 



SrO 

Ti02 

SrO 

Ti02 

FI 

-10.37 

13.47 

6.84 

3.10 

-3.53 

LM 

-0.44 

-4.93 

5.34 

-5.38 

4.90 

RI 

-10.81 

8.53 

12.18 

-2.28 

1.37 


3.3.3. Relaxed-ion slab calculations 


The results of the relaxed-ion slab calculations allow us to complete the picture of 
the fully relaxed flexovoltage response of a SrTiOa slab in the plate-bending limit. 
[The beam-bending case is easily recovered by multiplying the reported values by 
T = Cxx,xxl {Cxx,xx+Cxx,yy)- By using the calculated elastic constants of bulk SrTiOa, 
reported in Table we find r = 0.77.] A summary of the results is reported in Ta¬ 
ble The respective contributions of the bulk and surface are, overall, in line with 
the available order-of-magnitude estimates.^ The values shown in bold font, i.e., the 
total flexovoltage coefficients of the two types of slab, comprise the main result of this 
work. Note that they depart substantially from the corresponding bulk coefficient, 
confirming the dramatic impact of the surface structural and electronic properties on 
the electromechanical response of the system. In fact, the aforementioned response co¬ 
efficients are even opposite in sign depending on whether a SrO- and Ti02-terminated 
slab is considered. This is a remarkable result, as it means that an atomically thin 
surface termination layer can modify, and even reverse, the flexovoltage response of 
a macroscopically thick sample. This constitutes a rather drastic departure from the 
characteristics of other electromechanical phenomena (e.g., piezoelectricity), where 
the details of the surfaces typically become irrelevant in the thermodynamic limit. 

It is interesting to note that the surface shows an even larger termination de¬ 
pendence at the frozen-ion level, but with opposite sign. The LM contribution to 
j^surf jg large, and depends so strongly on the termination that its inclusion 

results in a voltage reversal, both in the Ti02- and SrO-type slabs. (By contrast, 
the LM contribution to the bulk flexovoltage coefficient is relatively minor, about one 
order of magnitude smaller than any other value reported in the table, and has little 
impact on the final results.) To illustrate the reason for such a strong dependence, 
a microscopic analysis of the surface relaxations is provided in Fig. 12 In the SrO 


case, the layer-by-layer decomposition of the induced dipole shown in Fig. 12 ’e) has 


an oscillatory behavior whose amplitude decays exponentially as a function of the 
distance from the surface; as a consequence, the surface layer clearly dominates the 
overall response.]^ Instead, for the Ti02-terminated slab shown in Fig. 12 f), the sur- 


from such a geometry do not necessarily reflect, therefore, the response of a realistic system; they 

are quoted here mostly for illustrative purposes. _ 

^ Interestingly, the structural relaxation pattern in the unperturbed state. Fig. |l2[ a), appears very 
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Fig. 12. Static and induced ionic relaxations at the SrTiOs surface, (a-b): Ionic relaxations in 
the unperturbed slabs (displacements from ideal bulk-like sites). Circles, squares, diamonds and 
triangles correspond, respectively, to Sr, Ti, 0(Ti) and 0(Sr) atoms. (Cations are indicated by 
empty symbols, oxygen atoms by filled ones.) Negative values indicate inward displacements (i.e., 
towards the slab center). (c-d): Displacements induced by a uniform strain of the type Syy — uExx ; for 
the two oxygen atoms in the Ti 02 layers, only one value (their average displacement) is shown, (e-f) 
Layer-by-layer decomposition of the lattice-mediated contribution to the induced surface potential 
offset. Vertical lines indicate the position of the SrO (solid) and Ti02 (dashed) atomic planes. 
(Adapted from Ref. |27| ) 


face layer responds with a positive dipole instead of a negative one, in sharp contrast 
to the “underdamped” oscillatory behavior in Fig. 12'e). This behavior is probably 


due to the alteration of the bonding network, which we speculate to be much more 
profound at the Ti02-type surface than at the SrO-type one, whereby the boundary 
atoms no longer behave as bulk-like but rather as a distinct chemical entity. 

Apart from the obvious relevance of the above observations to the physics of 
SrTiOa surfaces, the analysis of Fig. p^e-f) carries a general message that we have 
already anticipated in the above paragraphs. Any single atomic layer near the surface 


similar to the induced relaxation pattern under an applied tensile strain. This suggests that the 
former might be, in fact, rationalized as a response of the system to a large surface stress. 
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has a remarkably large contribution to sometimes of the same order as (or even 

larger than) the overall flexovoltage response of the slab. In fact, the total open- 
circuit voltage results from the subtle cancellation of many contributions of dissimilar 
physical nature. This implies that exceptional care is needed when dealing with 
flexoelectric phenomena, either when performing the calculations or when interpreting 
the experiments. 


4. Conclusions and outlook 


In this Chapter we have described the main advances in the first-principles theory of 
flexoelectricity that have taken place during the past five years. The progress that 
emerges from these pages is undoubtedly impressive - we are at the stage where the 
full flexoelectric response of real materials, including bulk and surface effects, can be 
calculated ah initio with great accuracy. Still, much remains to be done before the 
field can be regarded as mature. We discuss here several research avenues that we 
identify as being of pivotal importance for future progress. 


• Theory of the current-density response. The most fundamental and com¬ 
plete framework for the theory of flexoelectricity is the current-response formal¬ 
ism introduced in Sec. |2.2| Unlike the charge-response formalism summarized 
in Sec |2.5[ the current-response approach is capable in principle of resolving all 
independent components of the flexoelectric tensor. However, two issues remain 
to be settled in relation to this approach. First, direct methods for obtaining 
the current response functions of Eq. (131 by computing the linear response 
to a phonon of small but finite wavevector q have not yet been developed and 
tested. Once implemented, this would allow for a finite-difference calculation 
n(2.7''') ([l^^ and thence, the electronic contribution in Eq. (181. 


of the P, 


Ct,K,0 


Second, some aspects of the connection between the current-response theory and 
the theory of charge responses (including surface charges) remain to be clarified, 
as discussed in the context of Eq. ([50|) and following Eq. (fToTl). A solution of 


these two issues would help put the theory of flexoelectricty on a truly sound 
footing. 


• Analytic derivation of the q-expansions. The conceptual foundation of 
most of the material treated in this Chapter is a long-wave expansion of certain 
physical observables as a function of the wavevector q of an acoustic phonon. 
The calculations described in Sec.j^were performed by taking such a q-expansion 
numerically via finite differences, which is computationally cumbersome. Ide¬ 
ally, it would be best to perform the expansion analytically, i.e., to derive the 
DEPT equations that directly yield the wavefunction response to a strain gra¬ 
dient perturbation. This would also be desirable in the context of the direct 
current-density implementation sketched just above. When implemented in an 
existing DEPT code, such methods would allow for a more straightforward calcu¬ 
lation of flexoelectric properties of materials, and thus foster a more widespread 
application of these techniques within the research community. 
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• Application to complex materials. Our focus in this chapter has been on 
materials with cubic symmetry. Clearly a proper theory that also covers crys¬ 
tals of lower symmetry is strongly required. The extension of the theory to such 
materials will require attention not just to the proliferation of independent pa¬ 
rameters in the flexoelectric tensor, but also to subtle physical issues having to 
do, for example, with the anisotropic electronic screening that occurs when the 
symmetry is reduced. In the case of crystals that are piezoelectric (and possibly 
also polar), care will be needed to separate the higher-order flexoelectric from 
dominant piezoelectric (and possibly spontaneous) polarization response. The 
application to insulating ferromagnets or antiferromagnets should introduce no 
special difficulties in most cases, but may involve subtleties for magnetoelec¬ 
tric crystals or when spin-orbit coupling is strong. A first-principles theory of 
flexomagnetism has yet to be developed. 

• Compositional gradients. An electric polarization can also arise in the pres¬ 
ence of a compositional gradient, e.g., in Bai_a;Tia;03 films.l^ To our knowledge, 
a proper theory of such an effect is lacking. Since a compositional gradient gener¬ 
ally also entails a strain gradient, some care will be called for in separating these 
effects and computing them independently before combining the contributions 
to make physically meaningful predictions. 

• Connection to higher-level models. With the techniques described here, 
one can in principle calculate the fundamental flexoelectric properties of an ar¬ 
bitrary material. To use this information in real physical problems, however, 
one typically has to deal with many additional issues that are intractable by 
means of direct first-principles simulation: large samples with complex shapes, 
temperature effects, etc. It would be very desirable in this context to be able 
to extract the relevant physical parameters from the ab initio calculations, and 
incorporate them in some higher-level theory (e.g. atomistic, effective Hamil¬ 
tonian, or continuum) where length- and time-scale limitations are much less 
stringent. A successful attempt in this sense has already been reported!^ still, 
consistently incorporating the latest first-principles developments into macro¬ 
scopic theories remains an open challenge. For example, it would be of crucial 
importance, for a realistic description of the flexoelectric effect, to extract the 
relevant surface-specific properties from the density-functional calculations, and 
incorporate them into the higher-level model. Making progress in this direction 
will also promote a closer interaction between different communities working on 
flexoelectricity (continuum numerical modeling, Landau theory, etc.), which we 
believe would have a strong positive impact on the field. 

In summary, there has been dramatic progress in the development of a full first- 
principles theory of flexoelectricity. Several important challenges remain, as discussed 
above, but at least these have been identified, and solutions appear to be within 
reach. In any case, the development of the theory of flexoelectricity has already 
revealed many fascinating links to other, at first sight unrelated, research areas (e.g., 
the relationship to transformation optics, where the use of curvilinear coordinates 


54 


M. Stengel and D. Vanderbilt 


facilitates the solution of complex electrical engineering problems). We believe that 
more surprises are in store, and will progressively emerge while further progress is 
made along the above lines. As the study of flexoelectricity touches so many subhelds 
of condensed matter physics, we expect cross-cutting progress that will likely benefit 
the hrst-principles materials theory community at large. All in all, we look forward to 
the day when predictive calculations of flexoelectric responses can become a routine 
part of the tool-kit of hrst-principles computational materials theory. 
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